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ABSTRACT 


Earlier  experimental  work  has  been  extended  to  evaluate  the  effect  of 
moisture  on  the  Hugoniot  of  playa.  For  engineering  applications  the  Hugoniot 
of  moist  playa  can  be  predicted  with  sufficient  accuracy  from  the  Hugoniot  of 
dry  playa  and  water  and  the  assumption  of  pressure  equilibrium. 

Isentropic  release  data  were  obtained  for  moist  and  dry  playa.  The  steep 
release  curve  (in  the  P-V  plane)  from  high  pressure  implies  an  irreversible 
phase  change.  Some  low  pressure  data  in  the  elastic-plastic  region  are 
presented. 

A  theoretical  discussion  of  various  forms  of  the  Mie -Grime isen  equation 
and  the  physical  basis  of  asymptotic  statistical  models  is  presented. 

Shock  stability  is  reviewed.  Phase  transitions  in  which  AV  <  0  are 
classified  according  to  the  signs  of  the  slopes  of  the  coexistence  curves. 
Relative  slopes  of  Hugoniots  and  isentropes  in  the  mixed  phase  region  are 
calculated.  The  results  of  the  theoretical  discussion  are  applied  to  transitions 
in  bismuth,  iron,  and  quartz.  Agreement  of  values  of  dP/dT  deduced  from 
shock  data  and  measured  directly  are  good  for  bismuth  and  poor  for  quartz 
and  iron. 

Calculations  of  spherical  shock  propagation  in  a  hypothetical  medium  that 
undergoes  a  phase  change  are  presented.  The  calculations  show  qualitatively 
some  types  of  pulse  shapes  that  may  be  expected  in  a  transforming  medium. 

It  is  concluded  that  the  proper  treatment  of  phase  changes  is  an  outstanding 
problem  in  predicting  equations  of  state  for  earth  materials. 
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1-INTRODUCTION 


The  goal  of  this  program  is  to  measure  in  some  detail  the  equation  of 
state  of  Nevada  Test  Site  play  a,  to  extrapolate  the  results  to  pressure  and 
temperature  regimes  beyond  the  experimental  range  using  existing  theoretical 
methods,  and  to  examine  the  sensitivity  of  shock  propagation  in  spherical 
geometry  to  reasonable  variations  and  uncertainties  in  the  equation  of  state. 

hi  the  previous  year 's  work  Hugoniot  measurements  were  obtained  on  dry 
samples  of  playa  of  two  different  porosities,  of  initial densitities  1.55  and  1. 95 
g/cm  (crystal  density,  2.  66  g/cm  )  in  the  pressure  range  40  to  500  kbar. 

These  results  showed  several  interesting  features  :(1)  the  differences  in  thermal 
pressure  due  to  differences  in  initial  porosity  are  small  and  imply  a  small 
value  (  <  1)  of  the  effective  Gruneisen  parameter,  (2)  the  pressure-volume 
curve  appears  to  be  multivalued  in  the  200-  to  400-kbar  range.  This  result 
shows  that  a  simple  Mie-Griineisen  equation  of  state  with  F  a  function  of 
volume  only  is  inadequate  and  probably  implies  the  existence  of  polymorphism — 
presumably  the  quartz-stishovite  transition,  which  is  known  to  occur  in  this 
pressure  range. 

During  the  current  effort  these  results  were  extended  in  three  directions: 

1.  The  effects  of  moisture  were  examined  by  measuring  Hugoniot  states 
of  samples  of  the  same  initial  (dry)  density  as  before,  viz.  ,1.55  and  1.95 

3 

g/cm  ,  but  with  approximately  10  percent  by  weight  of  water  added.  In  addi¬ 
tion,  samples  containing  as  large  a  fraction  of  moisture  as  possible  consistent 
with  controllable  sample  quality  were  tested. 

2.  The  experiments  also  determined  several  points  on  the  release  isen- 
tropes  from  shocked  states.  For  earth  materials  particularly,  experimental 
determination  of  these  curves  is  just  as  important  as  determination  of 
Hugomots  because  of  the  possibility  of  irreversible  phase  transitions  and  be¬ 
cause  of  uncertainties  in  the  proper  theoretical  treatment  of  the  effects  of 
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moisture.  Both  of  these  problems  severely  complicate  the  derivation  of 
isentropes  from  Hugoniots  so  that  customary  procedures  used,  for  example, 
for  metals  and  simple  ionic  solids  are  of  questionable  validity.  Some  un¬ 
expected  peculiarities  were  in  fact  observed  at  higher  pressures. 

3.  Preliminary  experiments  were  also  performed  in  the  very  low  pres¬ 
sure  range  to  investigate  the  pressure  region  in  which  compaction  to  crystal 
density  occurs.  Unfortunately,  insufficient  effort  could  be  devoted  to  this 
problem  to  give  clearly  reliable  results.  Some  measurements  were  obtained, 
but  the  results  should  be  regarded  as  tentative . 

The  experimental  methods  and  results  are  described  in  Section  3. 

In  Section  4  a  general  discussion  of  various  approaches  to  predicting 
equations  of  state  is  given.  Also  in  that  section  is  a  comprehensive  review  of 
the  thermodynamics  of  the  Mie-Gruneisen  equation  of  state  and  a  discussion 
of  the  different  forms  in  which  it  is  used  by  various  authors. 

The  physical  bases  for  theories  of  high  pressure  asymptotic  forms  of 
equation  of  state  are  reviewed  in  the  appendix  in  elementary  form  to  assist 
the  nonspecialist  to  understand  the  assumptions  underlying  these  theories  and 
to  assist  him  in  appreciating  the  difficulties  in  assessing  their  validity. 

Clearly  one  of  the  most  difficult  and  potentially  important  problems  in 
constructing  an  equation  of  state  is  that  of  phase  changes,  including  poly¬ 
morphism,  melting,  and  vaporization.  The  existence  of  polymorphism,  the 
location  of  phase  boundaries,  and  the  relevant  thermodynamic  parameters 
describing  transitions  must  at  present  be  determined  experimentally,  and 
measurements  in  most  cases  are  lacking.  Moreover,  the  effects  of  phase 
changes  on  Hugoniot  and  isentropic  forms  of  equation  of  state  and  on  shock 
propagation  has  so  far  received  little  attention. 

The  work  reported  in  Section  5  is  a  fundamental  and  general  treatment 
of  the  thermodynamics  of  phase  changes  with  particular  reference  to  their 
effects  on  the  Hugoniot.  Application  of  this  theory  to  existing  shock  measure¬ 
ments  in  iron  and  quartz  shows  substantial  discrepancies  between  theory  and 
experiment — possibly  due  to  nonequilibrium  effects. 
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In  Section  6  the  results  of  spherical  shock  calculations  for  an  equation- 
of-state  model  containing  the  major  elements  of  a  phase  transition  are  pre¬ 
sented.  These  show  the  qualitative  shock  structure  to  be  expected  for  a 
reversible  phase  change.  A  summary  of  the  results  of  parameter  variation 
studies,  including  the  previous  year's  effort,  is  also  given  in  that  section. 


-  - 
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2-SUMMARY 


This  report  describes  an  extension  of  work  reported  previously  on  the 
equation  of  state  of  playa  from  Area  5  of  the  Nevada  Test  Site.  As  in  the 
previous  effort  the  work  was  concerned  with  (1)  experimental  determination 
of  the  equation  of  state,  (2)  theoretical  interpretation  of  the  experimental 
data,  and  extrapolation  by  semitheoretical  means  to  pressure  and  temperature 
regimes  beyond  the  experimental  range,  and  (3)  shock  calculations  to  indicate 
the  sensitivity  of  spherical  shock  propagation  to  reasonable  variations  in  the 
equation  of  state. 

The  experimental  work  was  extended  to  include  the  determination  of 
release  isentropes  from  shocked  states  in  both  dry  and  moist  playa.  These 
curves  appear  relatively  uncomplicated  where  the  peak  pressure  is  less  than 
100  kbar,  indicating  only  some  degree  of  compaction  by  the  fact  that  the  free- 
surface  velocity  is  less  than  twice  the  shock  particle  velocity. 

The  release  curves  from  a  shock  pressure  of  about  270  kbar,  however, 
shows  two  interesting  features.  For  dry  playa  the  initial  slope  (high  pressure) 
of  the  isentrope  is  very  steep  in  both  the  P-u  and  P-V  planes.  The  fact  that 
it  is  steeper  than  the  Hugoniot  in  the  P-V  plane  is  clear  evidence  of  some 
form  of  change  of  state  since  such  behavior  for  a  simple  fluid  would  violate 
the  shock  stability  condition.  The  most  reasonable  explanation,  consistent 
with  other  independent  observations,  is  that  the  quartz  component  of  the  playa 
is  converting  irreversibly  to  stishovite.  At  lower  pressures  the  isentrope 

from  270  kbar  becomes  shallow,  possibly  due  to  reconversion  of  stishovite  to 
quartz . 

The  isentropes  for  moist  playa  also  appear  to  be  uncomplicated  from 
shock  pressures  of  about  100  kbar,  For  the  higher  shock  pressures  the  free- 
surface  velocity  is  appreciably  higher  than  twice  the  shock  particle  velocity. 
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It  seems  likely  that  this  is  due  to  vaporization  of  the  water  on  release  of 
pressure,  and  resulting  rapid  expansion  of  steam  ahead  of  the  surface  of  the 
solid.  The  experiments  also  determined  the  Hugoniot  curves  for  moist  playa. 
These  results  show  that  the  effect  of  water  on  the  Hugoniot  can  be  predicted 
accurately  enough  for  engineering  applications  by  assuming  that  the  solid  and 
the  liquid  are  shocked  along  their  respective  Hugoniots  and  that  pressure 
equilibrium  is  obtained.  The  question  of  thermal  equilibrium  is  thereby  ignored, 
and  the  model  is  clearly  oversimplified.  Nevertheless,  it  is  found  empirically 
that  satisfactory  agreement  is  obtained. 

A  few  shock  experiments  at  very  low  pressures  were  performed  to  investi¬ 
gate  the  region  in  which  compaction  to  crystal  density  occurs.  The  results 
should  be  regarded  as  tentative,  but  indicate  a  precursor  wave  of  about  0.1 
kbar  amplitude  traveling  at  a  velocity  of  0.5  km/sec.  More  thorough  investi¬ 
gation  of  this  pressure  range  should  be  conducted  before  definite  conclusions  are 
drawn. 

The  theoretical  work  during  this  period  presents  a  general  review  of 
approaches  to  the  problem  of  predicting  an  equation  of  state.  It  also  presents 
a  thorough  treatment  of  the  thermodynamics  of  the  general  Mie-Gruneisen 
formulation  and  shows  the  differences  in  the  forms  used  by  different  authors. 
Finally,  an  elementary  description  is  given  of  the  assumptions  upon  which 
asymptotic  high  pressure  and  temperature  forms  are  based. 

Because  of  the  evidence  for  polymorphism  in  the  solid  constituents  and 
vaporization  of  moisture  in  the  playa,  and  because  these  effects  cannot  now 
be  easily  treated  theoretically,  a  considerable  effort  was  devoted  to  the  effects 
of  a  phase  change  on  both  the  equation  of  state  and  on  shock  propagation. 
Comparison  of  the  predicted  Hugoniots  in  the  coexistence  region  for  iron  and 
quartz  with  experimental  measurements  shows  substantial  discrepancies. 

These  may  be  due  to  incorrect  interpretation  of  the  data,  or  to  nonequilibrium 
effects.  The  results  for  bismuth  agree  reasonably  with  theory. 

Shock  propagation  calculations  were  extended  to  include  in  a  qualitative 
way  the  major  features  of  a  phase  transition,  and  typical  pressure  profiles 
and  decay  curves  are  shown. 
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In  general,  the  two-year  program  has  established  the  Hugoniot  equation 
of  state  from  40  to  500  kbar  including  the  effects  of  variable  porosity  and 
moisture  content.  Some  release  isentropes  were  measured  and  preliminary 
measurements  in  the  very  low  pressure  region  obtained.  A  relatively  simple 
theoretical  equation  of  state  was  developed  that,  in  the  absence  of  phase  changes, 
appears  adequate  for  playa  and  perhaps  other  earth  materials.  A  major 
remaining  difficulty,  however,  is  the  prediction  and  proper  treatment  of  phase 
changes;  progress  was  made  in  the  application  of  thermodynamics  to  this 
problem.  Shock  calculations  for  a  simple  energy  source  and  spherical 
geometry  showed  that  peak  pressures  as  a  function  of  radial  distance  are  not 
highly  dependent  on  uncertainties  or  variations  in  the  equation  of  state  and 
some  insight  into  reasons  for  this  insensitivity  was  gained.  Of  potential 
importance  to  interpretation  of  field  measurements  are  the  pulse  shapes 
associated  with  phase  changes  because  it  is  often  tacitly  assumed  that  the 
peak  pressure  is  closely  associated  with  the  first  shock  arrival. 

Possible  directions  for  extension  of  this  work  include: 

1 )  More  thorough  investigation  of  the  very  low  pressure  range  where 
the  material  is  not  completely  compacted. 

2)  Further  investigation  of  phase  changes,  due  to  polymorphism  and  to 
vaporization,  theoretically  and  experimentally. 

3)  Model  tests  in  which  shock  propagation  and  decay  can  be  compared 
with  predictions  based  on  the  equation  of  state  as  established  thus  far. 
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3 -EXPERIMENTS 


G.  D.  Anderson,  J.  T.  Rosenberg,  and  A. 


L.  Fahrenbruch 


A.  INTRODUCTION 

The  purpose  of  the  experimental  program  to  be  discussed  is  to  gather 
shock  wave  data  on  Nevada  Test  Site  playa  which  can  be  combined  with  existing 
theory  to  yield  an  equation  of  state  suitable  for  machine  flow  calculations.  The 
current  experimental  phase  is  a  continuation  of  a  program  which  was  begun  in 
mid  1963.  The  explosive  systems  and  streak  camera  techniques  used  in  ob¬ 
taining  Hugoniot  data  were  described  in  an  earlier  report1  which  summarized 
the  work  at  the  end  of  the  first  year.  During  the  first  year  the  effort  was  di¬ 
rected  toward  gathering  Hugoniot  data  on  dry  playa.  It  was  found  that  it  was 
necessary  to  reconstitute  samples  by  pressing  sifted  soil  in  order  to  obtain  sat¬ 
isfactory  streak  camera  records.  The  native  material  contained  inhomogenei¬ 
ties  in  density  which  were  large  on  the  scale  of  the  experiments.  These  inhomo¬ 
geneities  lead  to  irregular  or  "ragged"  shock  fronts  passing  through  the  samples 
which  destroyed  the  necessary  precision.  Samples  which  were  pressed  from 
soil  which  had  been  sifted  proved  to  be  quite  satisfactory.  X-ray  analysis  and 

streak  camera  records  both  indicated  a  uniform  density.  The  initial  densities 
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of  the  dry  material  were  1. 95  g/cm  and  1.  55  g/cm  .  The  fully  completed 
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crystal  density  was  measured  to  be  2.  65  g/cm  .  The  densities  of  the  pressed 
samples  studied  were  greater  than  the  native  dry  density  of  the  soil  (1. 39g/cm  ) 
as  it  was  found  to  be  necessary  to  press  to  higher  densities  in  order  to  obtain 
samples  which  were  mechanically  strong  enough  to  be  used  in  experiments. 

The  two  densities  were  chosen  so  as  to  generate  two  Hugoniots  for  the  purpose 
of  estimating  the  role  of  thermal  pressure.  No  large  difference  between 


Hugoniots  was  observed.  In  the  course  of  the  present  work  the  measurement 
of  the  Hugoniot  of  dry  material  has  been  repeated  with  good  agreement  with 
earlier  data. 

The  work  of  the  first  year,  which  has  just  been  briefly  summarized,  was 
expanded  in  three  directions  during  the  second  year.  Each  of  the  three  new 
phases  brought  new  problems  which  required  technique  development.  The 
three  phases  were  the  measurement  of  the  Hugoniot  of  moist  playa  and  evalua¬ 
tion  of  the  effect  of  water,  the  study  of  release  isentropes  including  free -surface 
velocity,  and  the  study  of  the  low  pressure  behavior  of  dry  playa  in  the  1  kbar 
region  where  compaction  may  not  be  complete  and  nonhydrostatic  behavior  is 
expected.  The  Hugoniot  measurements  of  moist  playa  presented  the  fewest 
experimental  problems  as  the  tests  relied  heavily  upon  techniques  developed 
during  the  first  year's  effort.  The  isentrope  measurements  and  the  low  pres¬ 
sure  studies  involved  new  types  of  experiments.  As  the  work  on  these  phases 
progressed  it  became  clear  that  extensive  studies  would  require  an  effort  much 
larger  than  the  present  one.  However,  significant  progress  was  made  toward 
the  perfection  of  the  new  techniques  and  some  preliminary  data  were  obtained. 

B.  SAMPLE  PREPARATION 

During  the  current  phase  of  this  program  an  effort  has  been  made  to  im¬ 
prove  the  existing  techniques  of  dry  sample  preparation  and  to  develop  a  method 
of  constructing  high  quality  samples  of  uniform  moisture  content  and  density. 

For  reasons  already  discussed,  the  playa  as  it  is  received  from  the  field  is 
unsuitable  for  small  scale  shock  wave  tests. 

The  initial  step  in  preparing  a  soil  stock  from  which  to  construct  test 
samples  is  to  shake  the  soil  through  a  series  of  sieves.  All  that  passes  through 
a  No.  50  sieve  (297  micron  openings)  is  recovered.  The  material  at  this 
stage  contains  5  to  6  percent  moisture  by  weight.  A  portion  of  this  soil  is  dried 
in  an  oven  at  105°C  to  be  used  at  a  later  time  as  a  diluent  to  a  high  moisture 
content  stock  in  the  preparation  of  specimens  of  various  intermediate  moisture 
contents. 
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(1 )  Control  of  Moisture  Content 


To  prepare  homogeneous  high  moisture  soil,  a  weighed  amount  of  the 
ambient  soil  is  placed  in  a  400-  to  600-ml  beaker,  leveled,  and  covered  with 
four  or  five  layers  of  filter  paper  cut  to  fit  snugly  in  the  beaker.  It  is 
important  that  the  height  to  diameter  ratio  of  the  soil  does  not  exceed  one.  The 
filter  papers  are  then  covered  with  several  layers  of  paper  towel.  Water  is 
added  to  the  paper  towel  which  absorbs  the  moisture  and  allows  it  to  slowly 
diffuse  down  through  the  filter  paper  into  the  soil.  No  more  than  8  to  10  g  of 
water  should  be  added  at  one  time  and  not  more  than  1/8  of  the  total  water 
should  be  added  in  any  12  to  24  hour  period.  After  addition  of  the  water  the 
beaker  is  sealed  and  allowed  to  stand  for  at  least  one  day.  Before  more 
water  is  added,  the  soil  is  poured  into  a  larger  beaker  and  thoroughly  stirred. 

It  is  then  replaced  in  the  smaller  beaker  and  covered  with  the  filter  paper  and 
towels  prior  to  adding  more  water.  This  slow  process  of  water  addition  is 
repeated  until  the  desired  moisture  content  is  achieved.  Upon  completion  of 
this  process  the  moist  soil  is  stored  in  stoppered  flasks.  Soil  samples  of 
intermediate  moisture  content  are  made  by  mixing  the  moist  material  just 
described  with  the  oven  dried  soil  in  the  appropriate  proportions  and  storing 
for  one  day  in  a  stoppered  flask. 

(2)  Sample  Pressing 

Since  the  low  density  samples  are  fragile  and  require  some  external  sup¬ 
port  after  removal  from  the  pressing  die,  they  are  pressed  in  steel  rings  of 
l/8-inch  wall  thickness.  Prior  to  pressing,  one  face  of  the  ring  is  covered 
with  0.0007-inch  Mylar  which  is  aluminized  on  one  side.  The  Mylar  is  bonded 
to  the  ring  to  form  a  seal  and  then  the  assembly  is  weighed.  The  Mylar  covered 
ring  is  bolted  in  the  pressing  die  so  that  the  ring  and  die  axes  are  parallel.  A 
predetermined  quantity  of  soil  is  then  poured  into  the  die  and  spread  uniformly 
with  a  leveling  tool.  At  the  time  of  pressing,  all  soil  contains  some  moisture 
since  it  has  been  found  that  samples  pressed  from  dry  material  crack  upon 
pressure  release.  Dry  samples  are  pressed  from  soil  initially  containing  a 
small  amount  of  moisture  and  then  dried  in  an  oven  at  105°  C  for  several  hours 
after  pressing.  The  density  is  controlled  by  pressing  a  weighed  amount  of  soil 
into  a  given  volume  fixed  by  a  series  of  stops  on  the  pressing  die. 


All  samples  are  weighed  immediately  after  pressing  and  those  to  be  di  ted 
are  then  placed  in  an  oven.  Moist  samples  are  sealed  in  the  rings  by 
0.0003-inch  Mylar  to  prevent  moisture  loss,  reweighed,  and  mounted  on  the 
driver  plate.  Moist  samples  are  not  made  until  just  prior  to  shooting  in  order 
to  minimize  moisture  loss  which  would  occur  during  long  periods  of  storage. 

Samples  of  density  less  than  1. 9  g/cm  are  pressed  in  steel  rings  which 
serve  a  dual  purpose.  They  provide  lateral  support  for  the  samples  which  at 
low  densities  are  relatively  fragile.  Also  they  are  a  convenient  support  point 
at  which  to  glue  the  aluminized  Mylar  which  covers  the  playa  to  hermetically 
seal  it  and  to  record  the  arrival  of  the  shock  wave.  The  pressing  process 
assures  that  the  playa  will  be  in  intimate  contact  with  this  aluminized  Mylar 
top.  Similarly,  care  is  taken  that  there  will  be  good  surface  to  surface  contact 
between  the  bottom  of  the  sample  and  the  lapped  2024  aluminum  driver  plate 
by  pressing  the  sample  to  a  thickness  several  mils  greater  than  that  of  the 
support  ring.  In  the  case  of  dry  samples,  which  may  be  kept  free  of  moisture 
during  delays  between  pressing  and  mounting  on  the  shot  assembly  by  storage 
either  in  an  oven  or  a  desiccator,  the  playa  could  be  mounted  in  direct  contact 
with  the  aluminum  driver  surface. 

g 

Samples  of  density  greater  than  1 . 9  g/cm  require  no  lateral  support  and 
are  pressed  free  standing.  In  the  shot  assembly  a  plastic  ring,  again  thinner 
than  the  sample,  is  used  as  an  anchor  to  which  to  attach  the  aluminized  Mylar. 
These  samples  are  not  stored  between  pressing  and  mounting,  hence  both  dry 
and  moist  specimens  are  attached  directly  to  the  driver  plate  surface  with  no 
intermediate  layer  of  Mylar. 
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C.  HIGH  PRESSURE  SHOCK  WAVE  EXPERIMENTS 


(1)  Hugoniot  Experiments 

The  Hugoniot  data  for  all  samples  were  obtained  by  the  impedance  match 
2 

method  which  is  described  quite  completely  in  Section  3.B  (2)  of  Reference  1. 
This  method  permits  determination  of  a  point  on  the  Hugoniot  of  an  unknown 
material  from  knowledge  of  the  shock  velocity  alone,  if  the  shock  is  introduced 
into  the  unknown  through  a  standard  material  whose  Hugoniot  and  relief  cross 
curves  are  well  known.  2024  Aluminum  was  used  as  a  standard;  its  Hugoniot 
and  calculated  isentropic  relief  curves  were  obtained  by  private  communication 
from  Dr.  R.  G.  McQueen  at  Los  Alamos  Scientific  Laboratory. 

The  2024  aluminum  driver  plate  used  as  a  standard  on  the  Hugoniot  experi¬ 
ments  is  nominally  8  inches  in  diameter  and  3/8  inches  thick,  the  dimensions 
varying  somewhat  from  shot  to  shot.  A  plane  shock  wave  is  induced  into  the 
driver  either  by  detonation  of  an  explosive  train  in  contact  or  by  impact  of  an 
explosively  driven  flying  plate  as  described  in  Section  3.B.  (1)  of  Reference  1. 
The  measurements  necessary  to  apply  the  impedance  mismatch  method  —  (free- 
surface  velocity  of  the  aluminum  driver  plate  and  shock  velocity  in  the  playa 
sample)  —  are  made  in  the  manner  described  in  Reference  1.  Some  detail  re¬ 
finements  have  been  made  in  order  to  attain  a  higher  degree  of  precision  in 
these  measurements.  For  example,  each  experiment  includes  two  independent 
measurements  of  the  free -surface  velocity  of  the  driver.  The  thickness  of 
shims  used  to  protect  gapped  mirrors  from  air  shock  has  been  reduced  from 
0.  006  to  0.  004  inches  with  the  result  that  pertubations  upon  the  values  of 
velocities  measured  by  these  mirrors  is  negligibly  small  after  corrections. 

No  samples  or  arrival  mirrors  are  located  at  the  center  of  the  aluminum 
driver  in  experiments  involving  flying  plates  since  it  has  been  noted  that  in 
some  cases  the  shock  wave  arriving  at  the  front  surface  of  the  driver  plate 
has  a  small  radially  symmetric  dimple,  in  both  pressure  and  shape.  The 
precision  of  the  playa  pressing  process  has  been  increased  during  the  course  of 
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the  project.  The  results  of  the  Hugoniot  measurements,  in  the  form  of  graphs 
and  tables,  are  presented  and  discussed  in  a  later  section  of  this  report. 

(2)  Release  Isentrope  Experiments 

The  problem  of  determining  release  isentropes  for  playa  necessitated  the 
development  of  new  techniques.  The  method  chosen  is  again  based  on  the  im¬ 
pedance  mismatch  principle.  A  shock  of  known  strength  in  a  standard  aluminum 
driver  plate  is  used  to  introduce  a  shock  into  the  playa  sample  in  the  same 
manner  as  in  the  Hugoniot  experiments.  However  in  the  adiabat  shots  a  ma¬ 
terial  of  lower  shock  impedance  than  soil  is  mounted  in  contact  with  the  front 
surface  of  the  soil.  As  porous  playa  is  of  relatively  low  shock  impedance,  the 
only  suitable  materials  of  lower  shock  impedance  are  liquids.  The  initial 
shock  propagates  through  the  driver  and  playa  as  before  until  it  reaches  the 
playa -liquid  interface.  There  a  rarefaction  is  reflected  back  into  the  soil,  and 
a  shock  is  transmitted  into  the  liquid.  If  the  Hugoniot  of  the  liquid  is  known, 
a  measurement  of  the  shock  velocity  is  sufficient  to  specify  the  state  in  the 
liquid  behind  the  shock.  As  the  boundary  conditions  require  continuity  of 
pressure  and  particle  velocity  at  the  playa -liquid  interface,  this  state  in  the 
liquid  must  be  a  pressure  and  particle  velocity  state  on  the  release  isentrope 
of  soil.  The  zero  pressure  point  on  the  release  isentrope  is  determined  by 
observing  the  free-surface  velocity  of  the  playa  which  is  constrained  to  remain 
at  essentially  zero  pressure. 

Experimentally  it  is  much  more  difficult  to  obtain  the  measurements 
necessary  to  determine  adiabats  than  Hugoniots.  All  adiabat  measurements 
are  made  after  the  shock  has  passed  through  the  playa  specimen  thereby 
increasing  the  number  of  uncertainties  which  may  enter  the  problem.  In  the 
Hugoniot  experiments  an  aluminum  driver  is  used,  whereas  in  the  adiabat 
shots  one  can  think  of  the  shock  being  introduced  into  the  liquid  by  a  playa 
driver. 

The  final  experimental  design  chosen  for  the  release  isentrope  experi¬ 
ments  is  shown  in  Fig.  3-1  and  will  be  described  below.  Each  assembly  yields 
a  Hugoniot  point  and  three  points  on  the  release  isentrope  from  that  Hugoniot 
point.  For  shots  in  which  the  explosive  is  in  contact  with  the  driver  plate,  the 
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FIG.  3-1  EXPERIMENTAL  ASSEMBLY  TO  MEASURE  RELEASE  ISENTROPES  OF  NTS  PLAYA 
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driver  thickness  is  3/4  inch.  For  shots  in  which  a  flyer  plate  with  a  free  run 
is  used  to  initiate  the  shock,  the  driver  thickness  is  1/8  inch.  The  reasons 
for  this  difference  in  driver  thickness  and  the  explosive  systems  used  in  each 
case  will  be  discussed  later  in  this  section.  The  smear  camera  trace  from 
slit  1  gives  two  independent  records  of  the  free -surface  velocity  in  the  driver 
at  approximately  the  same  radial  distance  from  the  center  of  the  driver  plate 
ar  the  observations  on  the  state  in  the  playa  and  liquids  are  taken.  For  each 
of  these  measurements  the  shock  arrival  at  the  free  surface  is  recorded,  and 
the  transit  time  of  the  free  surface  across  a  1 /8-inch  air  gap  is  measured.  The 
recording  surface  of  the  gapped1  mirror,  the  side  toward  the  driver  plate,  is 
protected  from  premature  arrivals  such  as  air  shock  by  a  0.004-inch  iron 
shim.  Since  the  Hugoniot  of  iron  is  known,  it  is  possible  through  application 
of  the  impedance  mismatch  method  in  an  iterative  manner,  to  correct  the 
observed  transit  time  for  the  presence  of  the  shim.  In  actual  calculations  the 
correction  is  small,  less  than  2  percent,  as  the  shim  thickness  is  small  com¬ 
pared  to  the  gap. 

The  points  defining  the  release  isentrope  are  taken  from  slits  2  and  3. 

The  playa  sample  of  diameter  2-5/8  inches  is  covered  with  aluminized  Mylar 
to  record  shock  arrival  and  planarity  at  the  front  surface  of  the  specimen.  On 
the  upper  two  quadrants  of  the  cell  are  1 /8-inch  gapped  free-surface  arrival 
mirrors,  protected  by  0.004-inch  iron  shims  as  before,  to  measure  the  playa 
free-surface  velocity.  The  lower  two  quadrants  of  the  cell  are  covered  by 
water  and  ethyl  ether,  both  transparent  liquids,  to  a  depth  of  1/8  inch.  The 
transparent  covers  of  the  liquid  cells  have  a  l/4-inch-wide  reflecting  stripe, 
protected  by  the  customary  0.  004-inch  iron  shim,  to  record  the  arrival  of  the 
shock  at  the  liquid  free  surface.  On  the  middle  slit  two  shock  arrival  mirrors 
are  mounted  on  tne  driver  surface  in  order  to  be  able  to  monitor  the  shock 
velocity  in  the  playa.  This  measurement  is  used  as  a  consistency  check  only 
since  knowledge  of  the  playa  Hugoniot,  determined  in  the  earlier  research 
period,  and  measurement  of  the  state  in  the  aluminum  driver  are  sufficient 
information  to  specify  the  state  in  the  playa  by  the  impedance  match  method. 
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The  gapped  mirrors  above  the  aluminum  driver  surface  are  supported  by 
means  of  1/8 -inch  hardened  steel  dowel  pins  whose  diameters  are  held  to 
tolerances  of  0.0001  inch.  The  dowel  pin -shim -mirror  assembly,  which  is 
glued  directly  to  the  driver  surface,  is  shown  in  Fig.  3-4  of  reference  1. 

The  final  uncertainty  in  gap  thickness  with  such  an  assembly  is  less  than 
0.0002  inch  and  hence  is  negligible. 

The  playa  specimens  are  supported  by  a  steel  ring  of  1  /8-inch  height  and 
wall  thickness.  Lucite  rings  of  the  same  diameter  and  wall  thickness  are 
mounted  concentrically  on  top  of  the  playa  support  ring.  The  Lucite  rings 
are  divided  into  quadrants  by  Lucite  cross  members.  This  assembly  is  hand 
lapped,  checked  for  parallelism  of  top  and  bottom,  and  held  to  maximum 
deviations  in  thickness  of  0.0005  inch.  The  entire  assembly  is  covered  on 
top  by  a  circular  piece  of  slide  glass  with  1/4-inch-wide  aluminized  stripes 
oriented  as  shown  in  Fig.  3-1  and  mounted  facing  the  playa  specimen.  Each 
of  the  two  liquid  cells  is  checked  for  leaks  between  quadrants  or  to  the  outside 
by  filling  with  air  at  a  pressure  of  at  least  10  psi.  Using  air  rather  than  liquids 
to  check  for  leaks  prevents  contamination  of  the  cells  and  reduces  the  possibility 
that  either  of  the  liquids  used  may  attack  any  component  of  the  cell  assembly. 

It  is  thought  that  water  may  cause  the  aluminized  Mylar  to  relax  or  deteriorate 
at  a  very  slow  rate,  and  similarly  that  ethyl  ether  may  attack  Lucite  at  an 
equally  slow  rate.  It  has  been  determined  that  neither  of  these  processes 
occurs  during  the  time  intervals  involved  in  the  course  of  firing  these  experi¬ 
ments. 

It  has  been  observed  that  on  some  previous  experiments  the  aluminized 
Mylar  covering  the  playa  has  pulled  away  from  the  playa  surface.  This  is 
thought  to  be  caused  by  the  fact  that  the  experiments  are  constructed  in 
temperature  controlled  environment  and  fired  at  the  test  site  where  the  ambient 
temperature  may  typically  be  20°F  higher.  Hence  the  gases  filling  the  pores 
of  the  sample  and  trapped  there  may  expand  and  lift  the  Mylar  from  the  surface 
of  the  playa.  To  avoid  such  situations  a  system  of  applying  an  overpressure 
of  approximately  one  pound  per  square  inch  to  the  top  surface  of  the  aluminized 
Mylar  has  been  developed.  This  is  accomplished  by  filling  the  two  liquid 
cells,  opening  passages  to  the  upper  two  quadrants,  and  applying  the  over¬ 
pressure  by  means  of  a  balloon  to  the  entire  Lucite  assembly  which  is  still 
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hermetically  sealed  from  the  outside.  It  is  apparent  from  visual  observations 
when  the  Mylar  has  been  pressed  into  contact  with  the  playa  surface,  and  hence 
it  is  possible  to  apply  only  the  minimum  overpressure  required.  This  over¬ 
pressure  is  in  all  cases  taken  to  be  so  small  as  not  to  effect  the  initial  densities 
of  either  liquid. 

A  typical  smear  camera  record  from  such  an  experiment  is  shown  in 
Fig.  3-2.  Fig.  3-2(a)  is  a  still  photograph  of  the  shot  face  with  the  image  of 
the  streak  camera  slits  exposed  over  it.  The  streak  camera  record  is  shown 
in  Fig.  3-2(b).  The  record  from  slit  1  yields  the  transit  time  of  the  aluminum 
driver  free  surface  across  a  1 /8-inch  gap.  The  record  from  slits  2  and  3 
yield  the  release  isentrope  data.  The  measurement  of  shock  velocity  through 
the  playa  from  slit  2  is  not  as  precise  as  the  shock  velocity  measurements 
from  shots  designed  to  determined  playa  Hugoniot  points.  This  loss  of  preci¬ 
sion  occurs  because  the  samples  used  for  isentrope  measurement  are  of  large 
diameter  to  permit  all  measurements  for  a  single  isentrope  to  be  taken  from 
the  same  sample  raised  to  a  uniform  Hugoniot  state  by  the  initial  shock.  Due 
to  the  large  diameter  of  the  samples,  shape  and  arrival  time  of  the  shock  at 
the  driver  playa  interface  on  slit  2  are  not  well  known.  However,  the  accuracy 
of  the  measurement  is  sufficient  for  a  consistency  check  on  the  state  in  the 
playa. 

The  measurement  of  playa  free -surface  velocity  is  complicated  by  the  fact 
that  playa,  being  a  porous  r  aterial,  is  subject  to  jetting  as  the  shock  arrives 
at  the  free  surface.  The  effect  of  such  jetting  is  to  cause  the  free -surface 
arrival  recording  mirrors  to  yield  jagged  free-surface  arrivals  with  a 
corresponding  high  degree  of  uncertainty  in  interpretation.  In  order  to  smooth 
the  jagged  arrivals  a  shim  is  mounted  directly  on  the  surface  of  the  playa. 

If  the  shock  impedance  of  the  shim  material  is  greater  than  that  of  playa  and 
if  the  release  isentrope  of  playa  from  a  doubly  shocked  state  is  not  significantly 
different  from  that  for  the  singly  shocked  playa,  then  it  can  be  shown  that  the 
shim  will  achieve  the  playa  free-surface  velocity,  through  a  series  of  wave 
reflections  at  the  playa -shim  and  shim-air  interfaces.  If  the  time  in  which  the 
shim  accelerates  to  the  playa  free-surface  velocity  through  the  wave  reflections 
is  small,  then  its  presence  will  have  a  negligible  effect  upon  the  value  measured 
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for  the  playa  free -surface  velocity.  On  slit  2  of  Fig.  3 -2(a),  a  0.003 -inch 
aluminum  shim  was  used  and  on  slit  3  no  shim  was  used  on  the  playa  free 
surface.  Aluminum  is  a  more  suitable  material  for  such  a  shim  than  iron  or 
steel  since  it  is  of  lower  impedance  thereby  introducing  a  smaller  pertubation 
on  the  state  of  the  playa.  Since  the  shock  and  rarefaction  velocities  of  alumi¬ 
num  are  higher  than  in  iron,  it  will  reach  equilibrium  in  shorter  time.  Also, 
fewer  wave  reflections  are  required  to  achieve  equilibrium  because  of  the 
closer  impedance  match  to  playa.  The  free-surface  record  from  slit  2  provides 
direct  experimental  examination  of  the  shim  effect.  The  shim  was  purposely 
cut  wider  than  the  free-surface  arrival  mirror  number  5.  The  shim,  labeled  7, 
can  be  seen  protruding  from  either  side  of  mirror  5  on  Fig.  3-2(a).  Along  the 
free-surface  quadrant  there  is  a  distinct  line  at  approximately  the  same  point 
in  time  at  which  the  free  surface  impacts  its  arrival  mirror.  This  is  due  to 
air  trapped  between  the  playa  free  surface  and  the  glass  cover  of  the  Lucite 
cell  luminescing  as  shocks  reflect  back  and  forth  raising  its  temperature. 

The  record  shows  quite  clearly  a  time  interval,  t  ,  between  the  luminescence 
due  to  the  arrival  of  the  accelerated  shim  and  due  to  the  unobstructed  free 
surface.  The  shim  arrives  later.  When  the  total  transit  time  of  the  shim  on 
the  free  surface  is  corrected  by  this  factor  there  is  very  satisfactory  agree¬ 
ment  between  the  two  free-surface  cells.  This  effect  will  be  discussed  further 
when  the  data  are  presented.  The  designations  A  and  B  are  used  on  mirrors 
number  4  in  Fig.  3 -2  (a)  to  point  out  that  the  aluminized  Mylar  on  the  playa  is 
being  observed  through  an  air  gap  and  a  liquid  cell  respectively. 

(3)  Explosive  Assemblies  for  Hugoniot  Release  Isentrope  Experiments 

The  explosive  assemblies  used  to  initiate  the  shock  in  the  driver  for  the 
Hugoniot  experiments  are  discussed  in  section  3.B.  (1)  of  reference  1.  The 
only  change  made  during  the  experimental  period  is  that  the  free  run  distance 
of  the  l/8-inch  stainless  steel  flying  plates  is  reduced  from  1-1/2  inches  to 
1  inch.  For  the  adiabat  shots,  however,  the  problem  of  attenuation  of  the 
shock  amplitude  with  distance  as  the  wave  progresses  through  the  experiment 
is  more  severe  since  measurements  are  made  over  twice  as  long  an  interval 
from  the  driver  surface  as  in  Hugoniot  experiments.  For  shots  in  which  the 
explosive  is  in  contact  with  the  back  surface  of  the  driver  plate  one  would  like 
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to  use  a  configuration  which  gives  a  long,  relatively  flat,  pressure  pulse  at  the 

3 

front  of  the  driver.  It  has  been  shown  in  this  laboratory  that  a  slowly  decaying 
pressure  pulse  can  be  obtained  in  an  aluminum  plate  from  an  explosive  train  of 
plane -wave  lens,  Comp  B,  and  Baratol,  if  the  ratio  of  the  thickness  of  Baratol 
to  that  of  Comp  B  is  2  to  1  and  the  aluminum  plate  is  at  least  1/2  inch  thick. 

Such  a  system  is  used  and  produces  a  shock  of  approximately  180  kbar  in  alumi¬ 
num  with  planarity  of  breakout  at  the  front  surface  on  the  order  of  0. 01  psec 
variation  along  a  3-inch  slit  length.  Shock  attenuation,  which  is  inevitable  due 
to  the  inherent  characteristics  of  an  explosive  detonation,  is  not  monitored 
directly,  but  by  Gregson's  report  should  not  be  a  source  of  difficulty. 

Wave  initiation  by  flying  plates  seems  desirable  for  two  reasons.  Higher 
shock  amplitudes  are  possible  than  for  in-contact  shots  since  the  flyer  plate 
receives  momentum  gradually  over  its  free  run  and  then  gives  it  up  rapidly  on 
impact  thereby  delivering  an  impulse  in  a  short  time  resulting  in  high  pressures. 
In  addition  to  higher  pressures,  it  is  in  principle  possible  to  achieve  flatter 
topped  pressure  profiles  via  the  impact  mechanism  of  a  flying  plate.  According 
to  hydrodynamic  theory  the  wave  should  be  perfectly  flat  until  the  trailing 
rarefaction  from  the  rear  of  the  flying  plate  overtakes  the  original  shock,  as 
discussed  in  reference  1 .  The  difficulty  in  designing  flying  plate  assemblies 
is  that  not  enough  information  is  known  to  accurately  compute  the  point  at  which 
the  overtaking  will  occur.  On  the  basis  of  early  results  in  the  adiabat  program 
it  is  felt  that  such  attenuation  was  taking  place  in  the  region  in  which  measure¬ 
ments  were  being  made.  Hence  a  new  system  designed  to  minimize  the 
possibility  of  attenuation  has  been  designed.  This  involves  increasing  the  ratio 
of  flyer  to  driver  plate  thickness  from  1/3  to  2  and  changing  the  flyer  material 
to  be  identical  with  that  of  the  driver.  Earlier  experiments  to  obtain  Hugoniot 
data  made  use  of  a  steel  flyer  with  an  aluminum  driver  plate.  Increasing  the 
flyer-to-driver-thickness  ratio  creates  two  problems.  As  the  driver  is  made 
thinner  it  becomes  more  difficult  to  machine  to  the  necessary  degree  of  flatness 
and  planarity,  and  as  the  flyer  is  made  thicker  it  becomes  more  massive  and 
hence  achieves  lower  velocities.  The  first  problem  is  really  one  of  economics 
and  has  been  met  simply  by  increasing  the  care  taken  in  the  machining  process 
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and  rejecting  all  unsatisfactory  plates.  The  second  problem  is  partially  solved 
by  changing  the  flyer  material  to  aluminum  which  is  approximately  1/3  the 
density  of  steel.  However,  since  the  shock  impedance  of  aluminum  is  less 
than  that  of  steel,  aluminum  must  be  accelerated  to  a  greater  velocity  than 
steel  to  produce  the  same  target  pressure  upon  impact.  The  maximum  pres¬ 
sure  attained  with  the  new  system  is  500  kbar  in  the  driver  whereas  700  kbar 
is  easily  attained  using  l/8-inch  steel  flyer  plates,  however,  500  kbar  is 
adequate  for  the  purposes  of  the  current  release  isentrope  program.  Previous 
work  has  shown  that  shock  attenuation  occurs  in  aluminum  flying  plate  experi¬ 
ments  earlier  than  is  predicted  on  the  basis  of  hydrodynamic  calculations. 

The  premature  attenuation  is  attributed  to  elastic  relaxation  due  to  elastic 
relief  waves  propagating  at  velocities  approximately  20  percent  higher  than 
hydrodynamic  values.  Taking  this  20  percent  velocity  increase  into  considera¬ 
tion,  time-distance  analysis  of  the  wave  propagation  in  flyer  and  driver  indi¬ 
cates  that  the  present  systems  should  be  free  of  attenuation  within  the  driver- 
playa-liquid  assembly.  The  fact  that  the  flyer  and  driver  are  of  the  same 
material  means  that  there  is  no  impedance  mismatch  at  this  interface  and 
hence  the  driver  may  be  made  as  thin  as  is  desired  without  new  disturbances 
originating  when  reflected  waves  from  the  playa-driver  interface  reach  driver- 
flyer  interface.  The  advantages  of  this  are  reflected  in  the  above  mentioned 
calculation  placing  the  attenuation  region  well  beyond  the  time  interval  in  which 
adiabat  measurements  are  made. 

Two  explosive  trains  are  used  with  the  new  flyer  plates.  For  intermediate 
pressures  the  Composition  B-Baratol  system  is  used,  and  for  high  pressures 
a  plane -wave  generator  and  an  HMX  pad  are  used.  With  the  first  system 
driver  pressures  of  265  kbar  are  reached  with  maximum  deviation  from  planar¬ 
ity  at  the  front  of  the  driver  being  0. 03  nsec  over  a  3  inch  diameter.  Because 
of  the  flatness  of  the  pulse  produced  by  this  particular  explosive  system  at  the 
back  of  the  flyer,  spalling  of  the  flyer,  which  could  introduce  premature 
attenuation,  is  unlikely.  The  HMX  system  gives  driver  pressures  of  500  kbar 
and  planarities  of  0.  01  fisec  over  a  3  inch  diameter  at  the  front  of  the  driver. 
The  observed  high  degree  of  planarity  is  very  satisfactory. 
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D.  DATA  AND  RESULTS 


(1)  Hugoniot  Experiments 

Hugoniots  for  NTS  playa  in  five  different  initial  states  have  been  experi- 

3 

mentally  determined.  Two  porosities  of  dry  playa,  P  =  1.55  g/cm  and 
3  ® 

Pq  =  1 . 95  g/cm  ,  were  selected  for  study  during  the  first  year  of  the  pro¬ 
gram  in  order  to  ascertain  the  effects  of  thermal  pressure  on  the  equation  of 
state.  Some  Hugoniot  work  for  dry  playa  of  these  densities  has  been  repeated 
in  this  experimental  period  to  serve  as  a  consistency  check.  As  playa  in  situ 
is  moist,  two  new  densities  of  playa  containing  approximately  10  percent  mois- 
ture,  pQ  =  1.71  g/cm  and  pQ  =  2. 14  g/cm  ,  were  studied.  These  two 
densities  were  obtained  by  requiring  that  the  samples,  in  addition  to  containing 
approximately  10  percent  moisture  by  weight,  have  the  same  pore  volume  as 

3 

the  corresponding  dry  samples.  Thus  if  a  moist  sample  of  density  2. 14  g/cm 

3 

were  to  be  dried,  the  resulting  sample  would  have  a  density  of  1. 95  g/cm  . 

3 

Finally,  playa  of  moisture  content  19  percent  and  initial  (wet)  density  1.55  g/cm 
was  examined  in  order  to  observe  the  effect  on  the  Hugoniot  of  having  the  pores 
filled  to  a  high  degree  with  water.  This  represents  the  highest  moisture  content 
it  was  possible  to  introduce  in  samples  of  wet  density  1.  55  g/cm  .  It  corre¬ 
sponds  to  56  percent  of  saturation.  It  should  be  remembered  when  comparing 
Hugoniots  of  dry  and  moist  samples  of  the  same  density  that  there  is  necessarily 
a  variation  in  pore  volume  between  the  two. 

The  results  of  the  Hugoniot  experiments  are  presented  in  Tables  3-1  to 
V  and,  as  pressure-particle  velocity  plots,  in  Figs.  3-3  to  7.  The  tables 
are  divided  into  direct -contact  and  flyer-plate  shots  on  the  basis  of  the  manner 
in  which  the  shock  is  introduced  into  the  driver.  It  should  be  noted  that  for 
some  of  the  very  low  pressure  shots  it  was  necessary  to  replace  the  standard 
aluminum  driver  plate  with  one  of  brass.  As  the  Hugoniot  of  brass  is  steeper 
than  that  of  aluminum,  application  of  the  impedance  match  method  will  show 
that  the  pressure  achieved  in  some  specimen  material  through  use  of  a  given 
explosive  material  and  a  brass  driver  is  lower  than  that  achieved  using  the 
same  explosive  material  and  the  standard  aluminum  driver.  Also  some  shots 
were  fired  in  vacuum  in  order  to  determine  if  air  in  the  pores  of  the  playa 
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Table  3-1 

IHGONIOT  DATA  R)R  NTS  PI.  AY  A 

=  1 .  33  i  0.0  2  g/cm 

Moisture  Content  (M.C. )  =  0  percent  +  0.2 


SHOT 

NO. 

DRIVER  DATA 

Pl.AYA  DATA 

Explosive 

System 

2024 

Aluminum 
Driver 
Pressure 
(  k  1)  a  r  ) 

Spec l me  n 
No. 

Shock 

Ve  1  oc  1 1  y 
(mm  ^  sec  ) 

Tree-Surface 
Velocity 
(mm /Msec  ) 

Particle 

Velocity 
( mm  Ul  sec  ) 

Pressure 
(  kb  a  r  ) 

K  i  ii a  1  V  <»  l  lime 
( i  m  ^  ■  £  ) 

Direct  contact 

1 . 1 

P-60  + 

1 

10, 584 

2 "  Comp  B 

28  5 

10 

4.01 

3.21 

2.08 

0.310 

A 15 

3.92 

-- 

2.09 

m 

0.300 

10,907 

2  Comp  B 

27  5 

5)6 

3.85 

3.  19 

2.05 

122 

0.  302 

10,605 

138 

32 

2.  58 

1.37 

1.23 

49 

0.338 

10,99(, 

-- 

137 

2. 49 

1.46 

1.23 

47 . 5 

0.327 

10, 606 

-- 

* 

16  f) 

mm 

1 . 9() 

0.852 

0.821 

24.  1 

0.367 

10,698 

-- 

* 

1M> 

mm 

1.86 

0.849 

0.819 

24 

0.359 

1  ''B  stainless  steel 

fiver  plates 

P-80  + 

10, 549 

3"  UNIX 

660 

12 

6.  38 

6. 13 

3.64 

357 

0.281 

11,053 

3"  I1MX 

602 

68 

6.35 

-- 

3.40 

331 

0.301 

11,131 

3"  KMX 

6 12 

71 

5.85 

.  „ 

3.  52 

318 

0.257 

69 

5.84 

-• 

3.  52 

318 

0.256 

11,173 

4"  UNIX5 

550 

7  t 

5.  54 

-- 

3.28 

282 

0.263 

1  1 ,  17  1 

Comp  B 

49  1 

75 

5.  38 

.  . 

3.02 

252 

0.283 

7  3 

5.34 

-  - 

3.04 

251 

0.279 

10,945 

2  Comp  B 

478 

36 

5.25 

4.79 

2. 99 

244 

0.277 

10,699* 

2"  Comp  B 

476 

8 

5.  14 

5.24 

2.97 

239 

0.271 

1 6 

5.  12 

-- 

2.98 

238 

0.270 

10, 58h 

2  Comp  B 

445 

A  9 

5.02 

4.90 

2.85 

220 

0.282 

* 

Brass  driver. 

^  Vacuum  shot. 

^  // 

>  16  stainless  steel  flyer  plate. 
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Table  3-11 


IIHGONIOT  DATA  FOR  NTS  PI, AY  A 

/>0  1.05  ±  0.02  g/rJ 

Moisture  (ion  ten  t  (M.C.  )  =  0  percent  +  0.2 


SHOT 

NO. 

DRIVKH  DATA 

PLAYA 

DATA 

Kx|»  l«a  i  ve 

System 

2024 

Aluminum 
Driver 
Pressure 
i k6a  r  ) 

Spe  c i me  n 
No. 

Sh  oc  k 

Ve  1  oc  i  t  y 
( mm  Msec  ) 

Free -Su  r fare 
Ve I oc  i  t  y 
( mm //Xscr  ) 

Particle 

Velocity 

Imm/fisec) 

Pressure 
( k  ba  r  ) 

Final  Vo  1 ume 

(c«Vg) 

l)i  rect  contact 

l'-I.O  + 

la  J-.  1 

10,467 

2"  Comp  li 

288 

E5 

4.55 

3.  13 

1.03 

172 

D1 

4.68 

-- 

1.92 

176 

10, 468 

158 

El 

3.35 

1 . 66 

1.26 

82 

0.  321 

D2 

3.34 

-- 

1.26 

82 

0.323 

10,  5  18 

16  3.4* 

D18 

2.57 

0.895 

0.772 

30 

0.355 

D12 

2.52 

0.775 

38 

0.354 

l^fl”  stainless  steel 

flyer  pi  ates 

P-80  + 

10, 585 

3"  IIMX 

610 

a3 

6.27 

4.76 

3.22 

394 

0.250 

10, 460 

2  Comp  l> 

50  3 

K7 

5.65 

-- 

2.84 

314 

0.  255 

10, 600 * 

2  Comp  H 

476 

4 

5.  18 

-- 

2.74 

205 

0.256 

*  Brass  driver. 
^  Vacuum  shot. 
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Table  3- III 

IIICONIOT  DATA  FOR  NTS  PI.  AY  A 
1.70  ±  0.01  K/Cm:< 

Moisture  Content  (M.C.  )  =  9.6  percent  t  0.4 


DRIVER  DATA 

PLAYA  DATA 

SHOT 

NO. 

Explosive 

Sy  s  tern 

2024 

Aluminum 

Driver 

Pressure 

Spec l me  n 
No. 

Shock 
Velocity 
( mm /Msec  ) 

Free -Surface 
Velocity 

( n’m/Msec  ) 

Par  t i c  1  e 
Ve 1 oc  1 1 y 

1  mm/isec  ) 

Pressure 
(  kb  a  r  ) 

Final  Vo  1  ii  me 

(cm  V  ) 

(  k  li  a  r  ) 

Direct  contact 

P-60  + 

Id,  784 

2  Comp  B 

285 

v  4 

4.  44 

3.39 

1.98 

1  50 

0.325 

i  1  ') 

4.  43 

-  * 

1.90 

149 

0.325 

10,605 

-- 

138 

14 

2.99 

1 . 56 

1.  17 

59.5 

0.358 

10,698 

-- 

166* 

3 

2.  36 

1.19 

0.795 

32 

0.  390 

10 , 606 

-- 

16  5* 

I 

2.  36 

0.950 

0.797 

32 

0.  389 

1/8  stainless  steel 

flyer  plates 

P-80  + 

10,549 

3"  UNIX 

660 

■) 

6.  58 

6. 46 

3.  50 

393 

0.274 

10,586 

2  Comp  B 

445 

<vf) 

5.  40 

-• 

2.7  2 

250 

0 .  29 1 

* 

lire  ss  dri  ver, 
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Table  3- IV 


IIUG0N10T  DATA  FOR  NTS  PI, AY  A 

f>Q  -  2.14  ±0.01  g  /  c  m  ^ 

Moisture  Content  (M.C. )  =  9.4  percent  i  0.2 


DHIVKH  DATA 

PLAY A  DATA 

SHOT 

NO. 

Kxp 1  os l v  e 

Sy stem 

2024 

Aluminum 
Dr i ve  r 
Pre  ssu  re 
( kb  a  r  ) 

Spec i me  n 
No. 

Shoe  k 

Ve 1 oc  i  t  y 
( mm/^sec ) 

Free-Surface 
Ve 1 oc  i  ty 
( mm/Msec  ) 

Particle 
Velocity 
(mm/Msec  ) 

Pressure 
(  kbar ) 

Final  Vo  1 ume 
(  cm  V g  ) 

Direct  contact 
l’-fiO  + 

mm 

10,  4(>7 

2  Comp  B 

288 

G8 

G3 

5.00 

4.72 

3. 18 

m 

196 

187 

119 

10, 4(>8 

-- 

1  S3 

G10 

G2 

3.80 

3.78 

1.67 

1.18 

1.18 

97 

05.5 

0.325 

0.321 

i  o ,  5  tfl 

-- 

If)  3* 

(i)  3 

3.52 

3.60 

1.36 

0.727 

0.724 

55 

55.9 

0.370 

0.372 

1/8"  stainless  steel 
flyer  plates 

P-80  + 

io,  r.8r) 

3"  MMX 

(ilO 

0)10 

6.81 

6.  26 

3.04 

444 

0.260 

10, 4d ^ 

2"  Comp  B 

303 

G6 

6.20 

5.01 

2.68 

359 

0.264 

* 

Brass  driver. 
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Table  3-V 

IIKiOMOT  DATA  FOR  NTS  I’l.AYA 
Pq  =  1 . 55  +  0. 0 1  g 'em  S 

Moisture  Content  (M.C.  )  =  18.9  percent  ±0.2 


DHIVKR  DATA 

PLAYA  DATA 

SHOT 

NO. 

K  x  p  1  o  s  l  v  e 

Sy  stem 

2024 

Aluminum 

Driver 

Pr  ensure 
(  k  b  a  r  ) 

Spec i me  n 
No . 

Shock 
Velocity 
(  mm  'Msec  ) 

Free -Sur  f  ace 
Velocity 
( mm/Msec  ) 

Particle 
Ve 1 ol  l  t  y 
(  mm//xsec  ) 

Pressure 
( kb  a  r  ) 

Final  Vo  1 u  me 
(  c  m  *  k  ) 

Direct  contact 

P-f>0  + 

10,00? 

2"  Comp  11 

27  5 

79A 

1.55 

2.53 

1 .  Of, 

138 

0  .  3(>(> 

1 0 , 8  <  1 0 

-- 

130 

45 

3. 22 

1.05 

1.17 

58.5 

0.410 

10, OOP 

-- 

137 

8  5 

3.00 

1.67 

1.  18 

57 

0.  300 

1  8"  stain  1  ess  steel 
flyer  plates 

P-80  + 

11,1 30 

3"  HMX 

(CIO 

70  A 

8  3A-  1 

().<)(> 

(..51 

-- 

3.  52 

3.  54 

3h() 

357 

0.  301 

0.  20  5 

11,0  33 

3"  UNIX 

(.02 

8  1  A-  1 

(,.08 

5.  Rh 

3.44 

324 

0.  280 

10,801 

2  Comp  11 

r  5 

32 

5.70 

(, .  50 

2.01 

257 

0.  315 

10,045 

2"  Comp  H 

478 

82 

5.  (.(. 

(,.20 

2.03 

257 

0.  311 
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PARTICLE  VELOCITY,  u  —  mm/jisec  M_S0M_I4 

FIG.  3-3  PRESSURE  vs.  PARTICLE  VELOCITY  IN  NTS  PLAYA 

(^q  =  1.55  g/cm3,  moisture  content  =  0  percent) 

affects  the  Hugoniot.  These  shots  are  appropriately  marked  in  the  tables.  On 
the  basis  of  the  results  it  appears  that  any  effects  due  to  air  in  the  pores  is 
smaller  than  experimental  error. 

The  densities  and  moisture  contents  which  specify  the  initial  conditions  of 
the  playa  are  recorded  at  the  top  of  each  of  the  tables.  The  quoted  tolerances 
in  densities  refer  to  the  maximum  deviations  that  were  actually  observed.  The 
value  of  the  average  density  of  any  given  sample  is  measured  to  within  1/4  per¬ 
cent.  The  tolerances  in  moisture  content  are  estimates  based  upon  observations 
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PARTICLE  VELOCITY,  u  —  mm//xsec  ob-So5»-i9 

FIG.  3-4  PRESSURE  vs.  PARTICLE  VELOCITY  IN  NTS  PLAYA 
{pQ  1.95  g  cm^,  moisture  content  0  percent) 

of  the  rates  at  which  dry  and  moist  playa  gain  and  lose  moisture,  and  upon 
known  variations  in  the  moisture  content  of  the  stocks  from  which  nominally 
similar  samples  were  pressed. 

Playa  free-surface  velocities  have  been  measured  in  many  of  the  Hugoniot 
experiments.  The  purpose  of  these  measurements  is  to  give  insight  into  the 
qualitative  behavior  of  the  free-surface  velocity  as  a  function  of  pressure. 
Inclined  mirrors  rather  than  gapped  mirrors  are  used  since  inclined  mirrors 
monitor  free-surface  velocity  continuously,  and  are  able  to  observe  any  struc¬ 
ture  which  the  shock  might  have  such  as  a  double  wave  induced  by  a  phase 
transition.  Iron  shims  are  used  on  the  free  surface  whereas  in  the  more 
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PARTICLE  VELOCITY,  u  - mm/usec 

r  06- SOW-  16 

FIG.  3-5  PRESSURE  vs.  PARTICLE  VELOCITY  IN  NTS  PLAYA 
(p q  -  1.70  g/cm^,  moisture  content  =  9.6  percent) 

sophisticated  adiabat  experiments  aluminum  shims  were  used  as  described  in 
a  previous  section.  Free-surface  velocities  measured  using  aluminum  shims 
indicate  that  the  free-surface  velocities  measured  earlier  using  iron  shims 
were  systematically  low.  The  variation  is  on  the  order  of  10  percent  at  100  kbar, 
and  less  than  2  percent  above  300  kbar.  This  variation  is  probably  due  to  shock 
attenuation  as  the  wave  passes  through  the  experiment  and  to  the  nonnegligible 
time  required  for  the  shim  to  accelerate  to  the  playa  free-surface  velocity.  In 
addition,  the  random  error  is  greater  for  measurements  made  with  inclined 
mirrors  than  gapped  mirrors  since  the  inclined  are  more  sensitive  to  shock  tilt 
and  curvature.  The  resulting  free-surface  velocity  measurements  are  useful  for 
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500 


FIG.  3-6  PRESSURE  vs.  PARTICLE  VELOCITY  IN  NTS  PLAYA 
(P  o  2.14  g/cm3,  moisture  content  9.4  percent) 

giving  the  general  nature  of  free -surface  velocity  dependence  on  parameters 
such  as  shock  velocity  or  pressure,  and  for  showing  the  relative  changes  in 
free-surface  velocity  with  the  variation  of  the  density  and  moisture  content. 

Figures  3-3  to  3-9  show  the  experimental  results  in  the  form  of  five  pres¬ 
sure,  particle-velocity  plots,  and  two  pressure  specific  volume  plots  for  playa 
in  the  various  initial  states.  Figure  3-3,  =  1.55  g/cm  ,  moisture  con¬ 
tent  =  0  percent,  has  error  brackets  on  two  typical  points.  These  are 

representative  of  the  errors  to  be  associated  with  points  in  the  high  and  low 
pressure  ranges  of  each  of  the  Hugoniots.  They  are  probable  random  error, 
rather  than  maximum  error,  calculated  from  estimates  of  the  uncertainty  in 
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Q^- - 1 - 1 - - - 1 

0  10  2  0  30  40 

PARTICLE  VELOCITY,  u  —  mm/usec 

GB-S059-IB 


FIG.  3-7  PRESSURE  vs.  PARTICLE  VELOCITY  IN  NTS  PLAYA 
(p  o  =  1.55  g/cm3,  moisture  content  -  18.9  percent) 

sample  density,  shock  path  length,  shock  curvature,  and  shock  transit  time 
as  recorded  by  the  smear  camera.  The  error  estimate  is  lower  than  in  the 
previous  year  partially  because  of  technique  refinements,  but  primarily  because 
of  the  availability  of  a  new  film  reader  having  considerably  better  resolution 
than  that  which  was  previously  attainable. 

The  five  Hugoniots  are  distinct  in  the  pressure,  particle^velocity  plane 
mostly  because  of  the  variation  in  initial  densities.  The  data  presented  in 
the  thermodynamic  plane,  pressure-specific  volume,  are  less  sensitive  to  the 
initial  density  and  the  various  Hugoniots  are  much  closer.  In  fact,  after  initial 
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400 


FIG.  3-8  HUG0NI0T  AND  RELEASE  ISENTROPES  FOR  NTS  PLAYA 
(l>  o  1.55  g  cm3,  moisture  content  -  0  percent) 


porosity  has  been  removed,  the  only  significant  difference  between  the  various 
samples  is  due  to  the  moisture  content  and  is  relatively  small.  It  is  interesting 
to  note  that  in  the  P-V  plane  it  is  possible  to  generate  the  Hugoniot  of  moist 
playa,  PQ  =  1.55  g/cm  ,  moisture  content  19  percent,  in  the  intermediate 
pressure  region  quite  closely  from  the  Hugoniots  of  water  and  of  dry  playa  of 
the  same  density.  This  is  done  by  making  the  simple  assumption  that  the  water 
and  playa  making  up  a  moist  sample  act  independently  and  are  at  the  same 
pressure  behind  the  shock.  Differences  in  temperature  are  neglected.  Since 
both  Hugoniots  are  known,  one  can  add  together  the  specific  volumes  of  playa 
and  water  at  various  pressures,  in  the  relative  proportions  which  each  are 
present  in  the  sample,  to  obtain  specific  volumes  of  moist  playa  at  those 
pressures.  In  this  way  a  Hugoniot  may  be  generated  which  agrees  quite  closely 
with  the  measured  one.  This  model  for  moist  playa  is  certainly  a  gross  over¬ 
simplification  of  the  actual  material,  hence  it  is  interesting  that  calculations 
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GB-5059-20 


FIG.  3-9  HUGONIOT  AND  RELEASE  ISENTROPES  FOR  NTS  PLAYA 
(Pq  1.55  g  "'em3,  moisture  content  =  18.9  percent) 


based  upon  it  are  not  greatly  different  from  actual  measurements.  A  similar 
result  was  reported  last  year  (1)  based  on  Sandia  data  for  dry  and  saturated 
sandstone. 


(2)  Adiabat  Experiments 

The  results  of  the  adiabat  measurements  are  presented  in  Table  3-VI  and 
Figs.  3-10  and  3-11.  The  table  is  again  divided  into  two  categories  according 
to  whether  the  wave  was  initiated  in  the  driver  by  a  flying  plate  or  by  explosive 
in  contact.  Release  isentropes  are  measured  for  only  two  of  the  five  initial 
playa  states  which  were  examined  during  the  Hugoniot  experiments.  It  is  felt 
that  the  other  initial  states  would  not  yield  significantly  different  results  and 
hence  did  not  warrant  the  additional  effort.  The  two  states  examined  are  both 

O 

of  density  1.55  g/cm  ,  one  dry,  the  other  of  moisture  content  19  percent.  The 
tolerances  on  the  soil  sample  parameters  are  as  quoted  in  the  previous  section 
on  Hugoniot  experiments. 
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FIG.  3-10  RELEASE  ISENTROPES  FOR  NTS  PLAYA  [pQ  -  1.55  g/cm3,  moisture  content  -  0  percent) 
Crosses  Denote  Error  Estimates 
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ISENTROPES  FOR  NTS  PLAYA  (pg  1.55g/cm^»  moisture  content  18.9  percent) 


Each  shot  yields  four  points  defining  one  release  isentrope  as  described  in 
an  earlier  section.  In  the  release  isentrope  experiments,  the  Hugoniot  state 
from  which  the  release  occurs  is  not  directly  measured.  Since  the  Hugoniot 
in  the  pressure,  particle -velocity  plane  is  quite  well  defined  from  all  previous 
data,  in  the  present  experiments  it  is  inferred  from  a  measurement  of  the  state  in 
the  aluminum  driver.  The  three  release  states  are  determined  by  measurement 
of  the  shock  velocity  through  the  liquid  reflectors  and  the  free-surface  velocity 
of  the  playa.  Typical  tolerances  are  shown  in  Fig.  3-10  for  points  on  the  high 
and  low  release  isentrope.  These  tolerances  are  calculated  by  a  similar  method 
to  that  used  in  the  Hugoniot  program.  Rectangular  tolerance  bars  are  not  used 
since  in  this  case  one  wants  to  know  the  variation  possible  in  a  curve  centered 
on  the  Hugoniot  rather  than  at  the  origin,  and  the  projection  of  rectangular  bars 
in  that  direction  does  not  cover  the  entire  error  region.  The  tolerances  shown 
for  the  free-surface  velocities  are  the  deviations  from  the  average  of  the  two 
values  measured  on  each  shot. 

It  should  be  noted  that  the  intermediate  pressure  release  isentropes  for 

both  moist  and  dry  playa  contain  fewer  points  than  the  other  isentropes.  For 

the  dry  playa  both  liquid  cells  were  filled  with  water  as  the  seal  between  them 

was  ruptured  after  the  leak  testing  procedure.  Hence  there  were  three  water 

cells  fired  with  the  intermediate  pressure  driver  system.  All  of  these  cells 

recorded  the  same  shock  velocity  within  experimental  error.  This  shock 

velocity  is  anomalously  low  in  the  sense  that  the  resulting  release  isentrope 

points  are  not  credible  on  physical  grounds.  On  a  pressure-particle  velocity 

plot  the  isentrope  points  fall  at  lower  pressure  and  lower  particle  velocity  than 

the  Hugoniot  points  from  which  they  originate.  This  phenomenon  could  be 

explained  by  postulating  that  there  exists  a  double  wave  in  water,  and  that  the 

velocity  which  is  being  measured  is  that  associated  with  the  first  wave.  In 

4 

support  of  this  explanation  it  should  be  noted  that  the  Russians  have  reported 
a  phase  transition  in  water  with  the  velocity  of  the  first  wave  being  within 
10  percent  of  the  velocity  we  have  observed.  These  three  points  are  recorded 
on  the  graph  but  are  not  taken  into  account  in  sketching  the  shapes  of  the 
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release  isentropes.  As  four  points  are  not  sufficient  to  determine  the  structure, 
if  any,  of  the  isentropes,  they  are  assumed  to  be  smooth  curves. 

An  interesting  difference  exists  between  the  free-surface  velocities  for  the 
two  moisture  content  playas.  For  dry  playa  the  free-surface  velocity  is  always 
less  than  twice  the  particle  velocity  behind  the  initial  shock.  For  moist  playa  of 
the  same  density  the  free-surface  velocity  is  greater  than  twice  the  particle 
velocity  except  for  the  lowest  pressure  in  which  case  it  is  almost  exactly  twice 
the  particle  velocity.  Free-surface  velocity  data  for  the  moist  playa  may  in¬ 
dicate  either  an  expansion  to  a  volume  considerably  greater  than  the  initial 
volume  or,  upon  release,  the  water  may  be  vaporized  and  thus  produce  a  free- 
surface  velocity  much  greater  than  that  of  dry  playa.  Some  of  the  early  shots 
using  inclined  mirrors  to  record  free-surface  velocities,  indicated  that  the 
mirror  was  sustaining  two  impulses  such  as  could  be  delivered  by  the  moist 
playa  if  water  vapor  and  then  playa  struck  it  successively.  As  the  scope  of 
the  project  did  not  permit  further  examination  of  this  hypothesis,  it  should 
be  borne  in  mind  that  the  free-surface  velocities  recorded  here  refer  to  the  first 
material  to  arrive  at  the  recording  mirror.  The  velocity  is  characteristic  of 
moist  playa  and  is  reproducible  to  the  accuracy  shown  by  the  tolerance  bars 
on  the  graph. 

Since  the  impedance  match  technique,  making  use  of  thecontinuity  of  pressure, 
and  particle  velocity  across  an  interface  between  two  media,  is  used  to  measure 
both  release  adiabats  and  Hugoniots,  the  data  are  naturally  obtained  in  the  form 
of  pressure,  particle-velocity  states.  The  conversion  of  a  Hugoniot  pressure- 
particle  velocity  state  to  a  pressure-volume  state  is  quite  readily  achieved 
through  the  application  of  the  Rankine -Hugoniot  jump  conditions.  Pressure- 
particle  velocity  adiabat  points  cannot  be  as  readily  converted  to  pressure- 
volume  points.  This  is  because  the  transition  from  the  initial  Hugoniot  state 
to  a  state  of  lower  pressure  and  higher  particle  velocity  and  volume  is 
achieved  by  a  continuous  process  through  a  rarefaction  wave  rather  than  an 
essentially  discontinuous  jump  as  in  a  shock.  Consequently,  the  jump  condi¬ 
tions  relating  the  two  states  in  the  case  of  a  shock  must  be  replaced  by  an 
integration  between  the  two  states  which  involves  all  intermediate  states. 
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Applying  the  equations  for  isentropic  flow  the  volume  at  some  state  in  a  rarefac¬ 
tion  wave  relieving  the  material  from  a  shocked  state  is  given  by 

u 

V  "  V1  '  /  (dP/du)g 
U1 

where  and  arc  the  volume  and  particle  velocity  state  behind  the  initial 

shock.  If  the  isentropic  pressure,  particle-velocity  curves  are  known,  the 

derivative  (dP/du)  can  be  calculated  and  the  integral  evaluated. 

s 

Smooth  curves  have  been  drawn  through  the  experimental  adiabats  in  the 
P-u  plane  in  order  to  map  them  into  the  P-V  plane.  If  slightly  different 
curves  were  drawn  through  the  data  different  P-V  adiabats  would  result  so 
that  the  curves  shown  in  Figs.  3-8  and  3-9  are  somewhat  arbitrary.  However, 
some  quantitative  conclusions  can  bo  drawn  from  the  general  shapes  of  the 
curves.  The  adiabats  coming  from  the  higher  pressure  Hugoniot  points  indicate 
the  adiabat  is  quite  steep  in  the  P-u  plane  compared  to  the  Ilugoniot.  This 
behavior  is  also  apparent  in  the  P-V  plane.  Such  a  phenomenon  has  been 

5 

observed  also  in  tuff.  This  small  increase  in  volume  or  particle  velocity 
with  decreasing  pressure  upon  release  can  be  explained  by  assuming  a  poly¬ 
morphic  phase  change  occurring  at  the  higher  pressures.  As  approximately 
50  percent  of  the  playa  is  silica  it  is  reasonable  to  suspect  a  transition  to 
stishovite.  The  adiabats  releasing  the  material  from  the  lower  pressure  states 
behave  in  a  more  normal  manner.  The  calculated  adiabat  for  the  moist  playa 

( p  ~  1.55,  moisture  1 8.  9  percent  by  weight)  releasing  from  the  highest 

°  3 

pressure  point  indicates  an  extremely  large  specific  volume,  1.65  cm  /g,  upon 

release  to  zero  pressure.  This  value  of  the  volume  is  again  dependent  upon 

the  assumed  curve  in  the  pressure,  particle-velocity  phase  for  the  isentrope. 

However,  the  high  free-surfacc  velocities  observed  for  this  moist  material  at 

high  pressures  imply  a  zero  pressure  volume  larger  than  the  initial  specific 
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volume.  The  effect  of  moisture  in  the  soil  appears  significant  in  releasing 
from  the  higher  pressures.  This  behavior  may  be  due  to  vaporization  of  the 
water  as  the  pressure  is  released. 

Values  of  T  ,  Grlineisen's  ratio,  have  been  estimated  from  the  slopes  of 
the  Hugoniots  and  isentropes  of  dry  playa  at  their  point  of  intersection.  The 
value  for  the  lower  pressure  isentrope  is  T  -  1.3.  At  the  higher  pressure 
point  r  -  -  16  .  This  anomalous  value  results  from  the  possible  phase 
change . 

E.  LOW  PRESSURE  SHOCK  WAVE  EXPERIMENTS 

The  extremely  low  stress  levels,  less  than  1  kbar,  may  well  be  the  most 
important  stress  region  for  study  from  the  point  of  view  of  application.  In  the 
case  of  a  blast  occurring  in  or  near  the  earth  the  majority  of  the  medium 
affected  by  the  ensuing  wave  motion  will  be  subjected  to  stresses  in  this  regime. 
Some  preliminary  dynamic  data  on  the  behavior  of  dry  playa  of  initial  density 
1.55  g/cm  were  obtained.  Since  the  techniques  of  inducing  very  low  ampli¬ 
tude  waves  into  soil  samples  and  recording  the  amplitudes  and  velocities  were 
unlike  any  techniques  used  in  other  phases  of  this  program,  the  largest  part  of 
the  effort  went  into  technique  development.  The  low  amplitude  waves  were 
induced  by  a  low  velocity  gas  gun  projectile  and  the  stress-time  recording 
was  done  with  a  quartz  pressure  transducer. 

(1)  The  Gas  Gun 

The  gun  consists  of  a  smooth  bore  2-l/2-inch-inside-diameter  barrel  which 
is  evacuated  ahead  of  the  projectile.  The  projectile  is  accelerated  down  the 
barrel  by  gas  introduced  from  a  high  pressure  reservoir.  Carefully  spaced 
electrical  pins  measure  the  projectile  velocity  near  the  target  assembly  as 
shown  in  Fig.  3-12.  The  mass  and  length  of  the  2-l/2-inch-diameter  pro¬ 
jectile  are  variable  and  may  bo  chosen  to  accommodate  each  experiment.  At 
the  present  lime  the  maximum  projectile  velocity  is  approximately  0.7  mm/ 
fj sec.  In  the  present  study  a  low  velocity  of  about  0.3  mm/ysec  or  100  ft/sec 
is  necessary.  The  main  problems  encountered  using  slow  projectiles  are 
velocity  control  and  tilt.  The  tilt,  a  measure  of  the  deviation  from  parallel 
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of  the  target  and  projectile  surfaces  at  impact,  is  measured  by  four  anival 
pins  in  the  target  assembly.  These  pins  are  shorted,  giving  rise  to  electrical 
pulses,  by  the  projectile  as  it  strikes  the  target.  The  pulses,  which  are 
binary  coded  to  assure  later  identification,  are  displayed  in  sequence  on  an 
oscilloscope.  The  time  interval  between  the  first  and  last  closure  measures 
the  total  tilt. 

(2)  The  Quartz  Pressure  Transducer 

The  quartz  gage  consists  of  an  x-cut  quartz  disk  which  has  a  conducting 
layer  evaporated  on  both  flat  surfaces.  A  circular  groove,  which  is  con¬ 
centric  with  the  disk  and  is  called  the  guard  ring,  is  machined  into  one  face 
and  divides  the  disk  electrically  into  two  regions.  Only  the  portion  of  the 
quartz  gage  within  the  guard  ring  is  used  for  recording.  The  outer  portion 
is  to  minimize  edge  effects  and  maintain  a  uniform  electric  field  in  the 
recording  area. 

When  a  pressure  pulse  traverses  the  quartz,  a  voltage  proportional  to  the 

difference  in  stresses  at  the  two  faces  is  generated.  The  gage  is  calibrated 

to  a  pressure  of  25  kbar  and  records  for  a  time  interval  equal  to  the  transit 

time  of  a  wave  through  the  crystal.  The  sensitivity  of  the  gage  is  about 

0. 8  volts/kbar.  A  more  comprehensive  treatment  of  the  behavior  of  the  quartz 

,  6 

gage  is  given  by  Graham  et  al. 

(3)  Shot  Assembly  and  Data 

3 

Dry  playa  samples  of  initial  density  1.55  g/cm  prepared  in  the  manner 
described  earlier  were  used  in  the  low  pressure  studies.  The  samples  were 
2-1/2  inches  in  diameter  and  approximately  1/8  inch  thick.  The  shot  assembly 
which  is  mounted  over  the  end  of  the  gun  barrel  is  shown  in  Fig.  3-13.  The 
projectile  strikes  an  aluminum  driver  plate  inducing  a  pressure  pulse  into  it. 
The  driver  plate  transmits  the  pulse  to  the  sample.  Upon  reaching  the  soil- 
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FIG.  3-13  DETAILS  OF  SHOT  ASSEMBLY  SHOWING  SAMPLE  AND  GAGE 


£ 


quartz  interface  the  pressure  pulse  is  reflected  back  into  the  soil  and  trans¬ 
mitted  to  the  quartz.  During  the  time  of  passage  of  the  initial  pulse  through 
the  quartz,  the  pressure  at  the  quartz -soil  interface  is  recorded  by  the 
transducer . 

An  aluminum  driver  plate  was  used  rather  than  impacting  the  projectile 
directly  into  the  sample.  Since  the  gun  barrel  is  evacuated  this  arrangement 
is  much  simpler  experimentally.  Very  low  projectile  velocities  were  used 
so  as  not  to  exceed  the  Hugoniot  elastic  limit  of  the  driver  plate.  If  a  double 
wave  system  emerges  from  the  aluminum  the  interpretation  of  the  gage 
record  becomes  more  confused  in  looking  for  a  double  wave  system  in  the 
soil.  The  aluminum  driver  is  equipped  with  pins  to  measure  projectile  tilt 
upon  impact.  A  layer  of  aluminum  foil  is  placed  on  top  of  the  soil  to  provide 
an  electrical  connection  to  ground  for  the  gage.  The  gage  is  placed  on  the 
soil  sample  and  the  entire  assembly  is  potted  in  C-7  epoxy  which  has  been 
doped  with  glass  beads  50  -  150^ in  diameter.  This  potting  is  done  to  more 
closely  match  the  impedance  of  the  surroundings  to  that  of  the  soil  and  gage, 
thereby  minimizing  edge  effects.  An  appropriate  resistor  is  added  from  tl^e 
guard  ring  to  ground  to  equalize  the  electric  field  in  the  quartz  betweerf  the 
guard  ring  and  center  electrode. 

Working  at  low  projectile  velocities  introduces  experimental  problems 
which  are  not  serious  at  high  velocities.  The  most  serious  of  these  problems 
is  tilt.  Since  the  projectile  velocity  is  much  lower  than  the  induced  pressure 
pulse  velocity,  large  refraction  effects  in  the  wave  front  arise  due  to  the 
nonsimultaneity  of  impact.  The  low  velocities  also  require  extreme  precision 
in  flatness  of  the  colliding  surfaces.  An  effort  was  made  to  maintain  all 
impacting  surfaces  flat  and  parallel  within  ±0.0001  inches.  It  was  found  that 
bowing  of  the  target  assembly  due  to  the  pressure  difference  when  the  gun 
barrel  is  evacuated  produces  a  nonplanarity  much  larger  than  the  tolerance 
specified  above.  To  eliminate  this  effect  an  auxiliary  vacuum  system  was 
added  to  the  back  of  the  target.  The  other  difficulty  in  working  at  low  projectile 
velocity  is  the  reproducibility  of  the  velocity  itself.  Frictional  forces  between 
the  barrel  and  projectile  become  quite  important  and  lead  to  wide  scatter  in 
velocity  for  the  same  initial  accelerating  gas  pressure  in  the  reservoir.  A 
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series  of  thirteen  shots  with  no  targets  was  fired  to  study  this  problem,  it  was 
found  that  using  a  massive  projectile  and  argon  rather  than  helium  as  an  accel¬ 
erating  gas  considerably  increased  the  reproducibility. 

The  pressure-time  profiles  recorded  by  the  quartz  gages  for  two  shots 
(11,467  and  11,468)  are  shown  in  Fig.  3-14.  Time  is  increasing  to  the  right. 
Both  records  were  subject  to  considerable  tilt  despite  the  precautions  taken  to 
minimize  it.  These  *>  shots  were  exploratory  and  were  the  only  experiments 
of  this  type  perform  Consequently,  in  the  absence  of  any  other  data  for  direct 
comparison,  these  results  must  be  regarded  as  tentative.  The  oscilloscope  rec¬ 
ord  for  Shot  11,467  shows  no  definite  double  wave  structure.  *  The  abrupt  change 
in  slope  in  the  rise  of  this  pulse  corresponds  to  a  pressure  of  about  0. 09  kbar  in 
the  quartz.  The  peak  is  at  about  0.  45  kbar  in  the  quartz.  If  an  elastic  precursor 
in  the  soil  is  present,  but  obscured  by  the  slow  rise  time  due  to  tilt,  it  would 
have  an  amplitude  of  0. 04  to  0. 09  kbar.  The  record  from  Shot  11,468,  given  by 
the  upper  trace,  clearly  shows  a  rise,  followed  by  a  plateau  and  then  a  second 
rise.  The  signal  on  the  lower  trace  is  from  the  guard  ring.  The  first  and  sec¬ 
ond  amplitudes  from  this  record  correspond  to  pressures  of  0. 12  and  0. 42  kbar 
in  quartz ,  respectively.  Interpreting  the  record  as  indicating  a  double  wave 
structure  in  the  soil,  the  first  wave  amplitude  would  be  0. 06  to  0. 12  kbar. 

Estimates  of  the  wave  velocities  may  be  obtained  from  the  time  interval 
between  impact  of  the  projectile  on  the  driver  and  the  arrival  of  the  wave  at 
the  soil -quartz  interface  as  indicated  on  the  gage  record.  The  transit  time  of 
the  input  wave  through  the  aluminum  driver  must  be  subtracted  out.  This 
transit  time  can  be  computed  from  the  known  Hugoniot  of  aluminum.  The  pres¬ 
sure  behind  the  second  wave  can  then  be  estimated,  ignoring  the  initial  wave 
which  is  small,  from  the  wave  velocity  and  the  fact  that  the  pressure,  particle- 
velocity  state  behind  the  second  wave  must  lie  on  a  relief  cross  curve  of 
aluminum.  The  pressures  and  particle  velocities  behind  the  main  wave  in  the 
soil  calculated  in  this  manner  are  presented  in  Table  3-VII.  The  pressures 
obtained  from  the  quartz  gages  for  the  main  wave,  which  must  be  the  pressure 
behind  the  reflected  shock  in  the  soil  at  the  quartz  interface,  are  lower  than 
would  be  expected.  As  the  impedance  of  x-cut  quartz  and  aluminum  are  quite 
close  at  these  low  pressures  it  would  seem  that  the  pressure  in  the  quartz 


♦The  amplitude  of  the  first  wave  was  taken  as  half  of  the  initial  pressure  rise. 
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would  be  at  least  double  the  initial  pressure  in  the  soil  if  the  soil  impedance 
remained  the  same  behind  the  shock.  If  the  main  wave  crushed  or  compacted 
the  soil,  one  would  expect  the  impedance  to  increase  and  the  pressure  to  more 
than  double  upon  reflection  from  quartz.  The  results  from  the  quartz  gage 
indicate  a  pressure  less  than  twice  the  estimated  initial  pressure,  implying 
that  the  reflected  shock  Hugoniot  of  the  soil  is  of  smaller  slope  than  the  initial 
Hugoniot  in  the  pressure,  particle -velocity  plane.  These  results  do  not  con¬ 
form  to  any  other  data  and  since  they  are  preliminary  they  must  be  regarded 
as  quite  tentative  pending  further  investigation. 


Table  3 - V I  I 

i.OW  PRESSURE  DATA  TOR  DRY  NTS  PLAY  A 


SHOT 

NO. 

SAMPLE 

DENSITY 

PROJECTILE 

VELOCITY 

(ft  sec  ) 

PROJECTILE 
TILT 
(y  set) 

WAVE  VELOCITIES 
(mm  usee  ) 

APPROXIMATE 
PRESSURE  IN  SOIL 
(  k  t>  a  r  ) 

APPROXIMATE 
PARTICLE 
VELOCITY  BEHIND 
SECOND  WAVE 
( mm/ Usee  ) 

<  £  «  m  S 

First 

Second 

First 

Second 

li,  it,: 

1 . 54 

If)  3 

2.6 

0.  37  4  0.04 

0.34  i  0.04 

0.  04-0.00 

0.30 

0.047 

1 1 ,  46R 

1 .  33 

143 

3.  1 

0.42  i  0.04 

0.30  ±  0.04 

0.06-0. 12 

0,23 

0.041 
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4 -THEORY 
Christian  Peltzer 


A.  INTRODUCTION 

When  attempting  to  obtain  equations  of  state  for  porous  earth  media  one 
faces  several  problems  in  addition  to  those  encountered  with  simple  solids, 
namely: 

a.  Already  in  its  initial  state,  the  medium  is  generally  a  three - 
component  system  consisting  of  a  solid  phase,  a  gas  phase 
(air)  and  a  liquid  phase  (water). 

b.  The  solid  phase  itself  is  an  essentially  isotropic  mixture  of 
many  different  compounds  and/or  mineralogic  species. 

c.  The  individual  constituents  of  the  solid  phase  are  often  com¬ 
plex  compounds  rather  than  simple  monatomic  crystallites. 

It  is  customary  to  treat  the  three  phases  as  noninteracting  systems  so  that  all 
thermodynamic  extensive  variables  are  obtained  additively  from  those  of  each 
phase  and  also  to  assume  that  pressure  equilibrium  is  attained  behind  the 
shock  front  and  maintained  during  pressure  release.  A  unique  specification  of 
the  thermodynamic  state  of  the  system  requires  of  course  further  assumptions, 
either  that  of  complete  thermodynamic  equilibrium  or  specific  assumptions  on 
the  behavior  of  each  phase  under  the  shock  transition  and  during  pressure 
release. 

Similarly,  the  solid  phase  is  often  treated  as  a  mixture  of  independent 
phases  and  its  thermodynamic  state  functions  are  then  calculated  additively 
from  those  of  the  simpler  individual  constituents,  usually  under  the  assumption 
of  complete  thermodynamic  equilibrium. 


Note:  See  List  of  Symbols  at  end  of  section, page  69. 
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A  meaningful  evaluation  of  the  validity  of  such  assumptions  would  require 
a  better  understanding  of  the  structure  of  a  shock  front  in  dense  mixtures  and 
little  can  be  added  presently  to  the  elementary  discussion  given  in  the  previous 
report. 1 

Hugoniots  for  wet  playas  obtained  as  described  there  from  the  Hugoniots  of 
solid  quartz  and  water  are  in  reasonable  agreement  with  the  experimental  wet 
playa  data.  Also  an  estimate  of  the  effect  of  a  decreased  Hugoniot  temperature 
in  wet  playa  (as  compared  to  dry  playa)  can  be  gained  by  taking  successively  for 
the  solid  phase  Hugoniot  the  experimental  data  of  solid  quartz  and  of  dry  porous 
playa. 

Calculated  Hugoniots  of  solid  mixtures  are  relatively  insensitive  to  the 

7 

averaging  procedure  used  to  obtain  them  from  the  individual  constituent 
Hugoniots  and  in  the  absence  of  a  better  understanding  of  the  physical  mechan¬ 
isms  underlying  the  propagation  of  shock  waves  in  solid  mixtures,  there  are 
little  real  grounds  for  selecting  any  one  of  the  several  averaging  schemes 

proposed  (per  weight  fraction,  per  molar  fraction, - )  in  preference  to  the 

others,  the  difference  between  the  resulting  Hugoniots  being  of  the  order  of 
the  experimental  uncertainty.  The  lack  of  sufficient  adiabatic  release  data 
precludes  at  the  present  any  analysis  of  the  same  problem  there. 

A  further  complicating  factor  in  analyzing  playa  data  is  the  existence  of 
polymorphic  phase  transitions  for  Si02,  in  particular  that  to  stishovite.  In 
view  of  the  particularly  large  volume  change  and  the  change  in  coordination 
number  associated  with  this  transition,  the  absence  of  conclusive  experimental 
evidence  for  or  against  its  occurrence  in  dynamic  compression  of  playa  makes 
any  comparison  of  the  experimental  data  with  postulated  equations  of  state 
rather  academic  for  the  time  being. 

Nevertheless,  one  may  conclude  from  a  comparison  of  the  experimental 
data  on  playas  and  quartz  and  from  the  limited  parameter  variations  performed 
with  the  equation  of  state  used  here,  that  porous  silica  is  a  satisfactory 
tentative  model  for  the  description  of  the  thermooynamic  behavior  of  the 
Nevada  playas  considered  in  this  program. 
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The  effects  of  small  variations  in  the  chemical  and  mineralogical  composi¬ 
tion  appear  of  little  significance  in  view  of  the  other  uncertainties  affecting 
both  the  experimental  and  theoretical  situations,  although  such  variations 
could  conceivably  be  of  more  importance  in  the  (little  studied)  low  pressure 
region  (<  20  kbars)  and  in  the  dynamics  of  eventual  phase  transitions. 

For  these  reasons,  the  theoretical  work  has  been  based  on  Si02  data 
only,  considered  as  a  single  thermodynamic  system  (see  reference  1). 

B.  DERIVATION  OF  EQUATIONS  OF  STATE  (EOS) 

Presently,  a  direct  first  principles  derivation  of  a  complete  EOS  for 
sufficiently  realistic  models  for  most  physical  systems  is  an  almost  impossible 
task  and  if  one  requires  an  EOS  to  be  valid  over  a  relatively  wide  range  of  the 
thermodynamic  variables  it  is  still  necessary  to  resort  to  semiempirical  for¬ 
mulations.  Most  equations  of  state  proposed  so  far  can  be  grouped  into  four 
types. 

a.  Purely  empirical  EOS;  these  are  obtained  simply  by  fitting 
largely  arbitrary  analytical  expressions  to  some  experimental 
data  such  as  a  Hugoniot;  if  such  EOS  may  be  of  interest  in  some 
calculations,  they  aro  but  an  economical  way  of  presenting 
specific  experiment?!  data;  they  are  totally  inadequate  for  any 
extrapolation  and  prediction  purposes  and  generally  of  no 
theoretical  value  or  particular  significance. 

b.  Semiempirical  EOS:  these  are  in  principle  obtained  by  choosing 
a  thermodynamically  consistent  functional  form  for  some  EOS 
and  by  determining  its  parameters  from  experimental  data  and/ 
or  theoretical  predictions  (see  c. ) 

c.  Theoretical  EOS:  These  are  obtained  from  first  principles 
according  to  the  laws  of  statistical  mechanics;  they  are 
presently  limited  to  a  few  simple  systems  (classical  and 
quantum)  such  as  a  free  electron  gas,  a  harmonic  lattice  solid, 
etc. ,  and  are  exact  within  the  well  defined  limits  of  the  model 
and  its  underlying  physical  laws. 
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d.  Semitheoretical  EOS:  In  this  disparate  group,  one  can  include 
many  EOS  which  although  partly  derived  from  first  principles, 
can  only  be  obtained  by  additional  assumptions  of  an  heuristic 
nature  whose  degree  of  validity  is  difficult  to  assess;  e.  g. , 
Thomas-Fermi,  (TF)  Thomas -Fermi-Dirac  (TFD),  etc. 

In  view  of  the  many  loose  statements  implying  the  contrary,  it  may  not  be 
superfluous  to  repeat  that  semitheoretical  EOS  such  as  the  T-F,  T-F-D,  etc. , 
are  in  no  way  exact  EOS  of  a  well  defined  model  and  that  at  the  present  there 
exists  no  true  estimate,  in  a  proper  mathematical  sense,  of  their  degree  of 
validity  in  any  particular  range  of  the  thermodynamic  variables,  but  only 
more  or  less  optimistic  guesses  as  to  their  applicability  (see  later). 

We  shall  limit  ourselves  here  to  a  consideration  of  some  semiempirical 
EOS  and  we  shall  also  discuss  a  few  questions  concerning  the  use  of  the  T-F, 
T-F-D  EOS  and  the  nature  of  certain  corrections  proposed  by  various  authors 
(see  also  Appendix. ) 

C.  SOME  THERMODYNAMICS  OF  EOS^ 

For  the  reasons  given  earlier  (see  also  Ref.  1),  we  consider  here  only 
single  thermodynamical  systems  (homogeneous,  isotropic)  describable  in 
terms  of  two  independent  variables,  taken  mostly  to  be  V  and  T  or  E*  and 
T,  where  E*  is  the  thermal  energy. 

(1)  The  General  Mie-Griineisen  EOS 

The  general  Mie-GrUneisen  EOS  is  defined  by  the  relation 

r=T(V,T)  (4-1) 

T  being  the  Mie-GrUneisen  ratio  defined  as 


^  For  notation  see  reference  1. 


(4-2) 
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This  EOS  was  discussed  in  reference  1  and  we  only  recall  its  main  features: 


(i)  it  is  an  incomplete  EOS  that  leaves  the  cold  energy  E„(V) 
unspecified 

(ii)  it  uniquely  det  rmines  all  other  adiabats  of  the  system  in 
the  V-T  plane,  these  are  the  characteristics  of  the  partial 
differential  equations, 

xr(V.  T>(|f  )v  -  v(H)  T  =  0  (4-3) 

S(V,  0)  =  constant  (=  0)  (4-3a) 


and  so  are  given  by 


z  =  0r(V,  T)  =  constant  (4-4) 

where  </>  j-(V,  T)  is  a  function  of  V,  T  uniquely  determined  by  T  (V,  T) 

(iii)  The  entropy  is  a  function  of  z  only,  arbitrary  up  to  general 
thermodynamic  restrictions  on  admissible  S(V,  T)  's  (positive 
definite,  etc. ) 

S  =  o(z)  (4-4») 

(iv)  A  complete  specification  of  the  system  requires  besides 

r (V,  T)  the  knowledge  of  E  (V)  and  that  of  o(z)  or  its 

c 

equivalent,  e.  g. ,  one  isobar  VpQ  =  VpQ(T)  or  one 

isotherm  PT  =  PT  (V)  (Trt  ^  0)  or  a  specific  heat 
*o  Ao  o 

C  v  =  Cy  (T)  ,  etc. 
vo  o 

Two  useful  equivalent  forms  of  the  general  Mie-Griineisen  EOS  are: 

P  -  PC(V)  -  T(lf-\,  -  (w)T  (4-5) 


55 


dEc 

dV 


(4-5a) 


P 


Ec  -  Ec(V,0) 


and 


P(V,T) 


r(v,T) 

V 


E(V,  T)  +  R(V,  T) 


(4-6) 


Rc<V> 


(4-6a) 


It  is  important  to  note  that  R  is  a  function  of  V  alone  if  and  only  if  T  does 
not  depend  explicitly  on  T  . 

(2)  Special  Cases  of  the  General  Mie-Grimeisen  EOS 

Two  special  cases  of  the  general  Mie-Griineisen  EOS  have  been 
frequently  used,  namely 

(i)  The  usual  Mie-Griineisen  EOS 


r  =  T  (V) 


usually  written  in  one  of  the  two  equivalent  forms: 


p  -  Pc(V)  =  ip  |e 


p 

C 


E  =  E  (V) 
c  c 


(4-7) 


(4-8) 


(4-8a) 


or 


+  R(V) 


p 


(4-9) 


(4-9a) 


4  = 


The  explicit  form  of  the  principal  thermodynamic  state  functions  for  this  case 
were  given  in  reference  1,  as  well  as  further  details  on  its  use  for  Hugoniot  and 
adiabat  calculations. 

(ii)  The  Hildebrandt  EOS 

In  any  domain  (T  ^  0,  T)^  v/here  E*  is  a  function  of  T  only,  the  general 
Mie-Griineisen  EOS  can  be  put  in  the  form 


P  -  PC(V) 


r<V,  T)  dE 


P 

c 


dE 

c 

dV 


E  =  E  (V)  E* 

c  c' 


T  ^  0 
E*(T) 


(4-10) 


or  equivalently 


PV  + 


dEc(V) 
d  In  V 


*  * 

E  =  E  (T) 


T  ^  0  (4-10') 


These  two  special  cases  are  not  mutually  exclusive  and  it  can  be  shown  that: 


The  usual  Mie-Griineisen  EOS  and  the  Hildebrandt  one  are  equivalent  in  a 
domain  (T  t  0,  T)  if  and  only  if  Cy  =  const  and  E*  is  of  the  form 
E*  =  CyT  +  const  ,  T  being  then  necessarily  a  function  of  V  alone,  given 
by 


r 


VS'T  (V) 
o 


(4-11) 


where  Stq(V)  is  the  entropy  along  the  Tq  -  isotherm, 
t Since  this  special  case  is  incompatible  with  the  3^  law,  Tq  must  be  >  0  . 
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x-l  . 


In  particular,  one  can  note  that  a  perfect  gas  and  a  Dulong-Petit  solid  both 
satisfy  this  equivalence  criterion. 

(iii)  More  generally,  if  F(V,  T)  is  of  the  form  Tc(V)rt(T),  explicit 
forms  for  all  thermodynamic  state  functions  can  be  obtained 
just  as  for  a  usual  Mie-Griineisen  system  Y  =  T  (V)  .  For 
other  forms  of  T  ,  the  characteristic  equation  of  (4-3)  cannot, 
in  general,  be  integrated  in  closed  form  and  particular  attention 
should  be  given  to  the  consistency  of  separate  assumptions  on 
the  forms  of  F  (V,  T)  and  other  thermal  entities  such  as 
E*  or  Cy  . 

(3)  Modified  Mie-Gruneisen  EOS 

Several  other  related  EOS's  have  been  considered  in  the  literature,  in 
particular  the  following  ones. 

a.  As  already  remarked,  an  EOS  of  the  form 


P(V,  T)  =  e(V,T)  +  X(V) 


(4-12) 


cannot  be  a  Mie-Griineisen  EOS,  i.e.,  the  function  f  (V,T)  cannot — if  it 
depends  explicitly  on  T—  be  identical  to  the  Mie-Gruneisen  ratio  T  defined 
by  (4-2).  The  two  functions  r(V,T)  and  F  (V,  T)  are  aciually  related  by: 


r-  r|  (ff)  -  v(|f) 


(4-13) 


It  is  possible  to  choose  T  such  that  Tq(V)  e  T(V)  but  then  one  must  also 
have 


R JV)  =  X(V)  =  P  (V)  - 


Tc(V)Ec(V) 


(4-14) 


58 


i.e.,  only  two  of  the  three  functions  X(V)  (or  R  (V)),  E  (V)  and  P  (V) 

/  \  \  c  /  c  c 

(or  rc(V)  can  be  chosen  arbitrarily,  the  third  one  being  given  by  Eq.  (4-14) 

A  detailed  study  of  this  modified  Mie-Griineisen  EOS,  (4-12),  can  easily 
be  carried  out,  but  will  be  omitted  here  since  this  form  does  not  appear  to 
offer  any  particular  advantages  over  the  general  Mie-Griineisen  EOS. 

8 

b.  Some  authors  have  made  use  of  the  following  EOS: 


P(V,  E)  -  E  +  X(V) 


(4-15) 


The  function  G(V,  E)  is  also  distinct  from  the  Mie-Griineisen  ratio  T  if  it 
depends  explicitly  on  E  ,  the  two  being  related  by 


r  -  g  = 


(4-16) 


If  (4-15)  is  to  hold  along  the  cold  isotherm,  one  obtains  the  equation 


dE 
_ c 

dV 


Ec(V)  +  X(V) 


(4-17) 


relating  G(V,  E)  ,  Ec(V)  ,  and  X(V)  . 

Q 

If  the  forms  of  G(V,E)  and  X(V)  have  been  chosen  a  priori  (with  adjust¬ 
able  parameters  determined  empirically),  Eq.  (4-17)  must  be  used  to  obtain 
the  cold  energy  Ec(V)  .  Such  a  procedure  is  perfectly  legitimate  for  semi- 
empirical  EOS  purposes,  although  it  may  be  preferable  to  choose  first  the 

forms  of  E  (V)  and  G(V,  E)  and  then  use  Eq.  (4-17)  to  determine  X(V)  . 
c 

E  (V)  is  a  relatively  accessible  quantity  both  theoretically  and  experimentally 
c 

and  this  procedure  would  make  easier  a  comparison  with  other  formulations 
and  a  verification  of  the  necessary  stability  requirements  for  E  (V)  .  The 

V  ^ 

more  serious  drawback  of  this  formulation  is  that  the  dependence  of  G  on 
Ec  and  E*  (and  so  on  T)  is  the  same  and  already  determined  by  the  relation 
(4-17)  along  the  cold  isotherm.  This  rather  artificial  feature  renders  any 
interpretation  of  the  thermal  dependence  of  G  difficult  and  precludes  any 
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nonpurely  numerical  comparison  of  this  EOS  with  other  formulation  and 
models. 

c.  For  these  reasons,  the  following  EOS  appears  preferable: 


P(V,E*)  -  Pc(V) 


_  G(V,  E*) 

V  ** 


P 

c 


dE 
_ c 

dV 


E 

c 


E  (V) 
c 


(4-18) 


The  ratio  G(V,  E*)  is  again  identical  to  the  Mie-Griineisen  ratio  T  if  and 
only  if  it  does  not  depend  explicitly  on  E*  ,  in  which  case  Eq.  (4-18)  reduces 
to  the  usual  Mie-Griineisen  EOS. 

In  contrast  to  Eq.  (4-15),  it  does  not  impose  any  a  priori  coupling  between 
the  cold  and  thermal  components  of  the  thermodynamic  -state  functions  and 
although  it  does  not  appear  to  offer  any  advantages  theoretically  over  the  general 
Mie-Griineisen  EOS,  it  can  be  of  interest  in  semiempirical  EOS  work  and  shock 
calculations. 

The  use  of  so  many  similar  but  distinct  EOS  (and  a  host  of  other  semi¬ 
empirical  and  semitheoretieal  ones)  is  unfortunate  inasmuch  as  it  makes  a 
comparison  of  various  authors'  results  a  tedious  and  difficult  process  while 
contributing  very  little  but  added  confusion  in  understanding  and  solving  the 
many  remaining  problems  in  this  area  The  general  Mie-Griineisen  EOS  (or 
eventually  the  modified  Mie-Griineisen  EOS  (4-18))  could  provide  a  sufficiently 
general  formulation  for  the  present  needs.  A  systematic,  unified  presenta¬ 
tion  of  the  existing  experimental  and  theoretical  data  for  the  many  systems 
invc  ited  up  to  now  would  provide  one  with  a  much  better  picture  of  the 
present  situation  and  help  considerably  in  identifying  those  areas  where  further 
work  is  most  needed. 


D.  DERIVATION  OF  A  SEMIEMPIRICAL  EOS 

When  seeking  a  semi-empirical  EOS  for  a  specific  system,  it  appears 
best  to  proceed  in  successive  steps  as  follows: 
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a.  First  choose  a  suitable  form  for  the  cold  energy  E  (V) 

c 

b.  Select  then  an  appropriate  thermodynamic  type  of  EOS,  e.g. , 
the  general  Mie-Griineisen  EOS 

c.  Starting  from  the  simplest  ones,  make  further  consistent 
assumptions  on  the  specific  form  of  the  thermodynamic  state 
functions  determining  the  thermal  components  of  the  system, 
e.g.,  T(V,T)  and  Cy(T)  . 

With  guidance  from  both  experimental  data  (static  compression,  Hugoniots, 
etc. )  and  theoretical  considerations,  considerable  progress  can  be  achieved  in 
an  orderly  and  systematic  way  towards  the  understanding  and  semiempirical 
representation  of  the  principal  factors  determining  the  EOS  of  many  materials 
over  a  fairly  wide  range  of  the  variables.  This  is  well  illustrated,  e.g.  ,  by 
the  work  of  the  Russian  school  (Kormer,  Al’tshuler  et  al. )  on  metals  (solid  and 
porous),  recently  extended  to  the  alkali  halides,  which,  if  still  unsatisfactory 
in  several  aspects,  appears  nevertheless  to  have  led  to  the  most  successful 
EOS  for  those  materials  presently  available. 

Most  purely  theoretical  results  are  limited  so  far  to  a  few  very  simple 
systems  which  provide  in  this  context  mainly  very  high  density  and  or  very 
high  temperature  asymptotic  models.  The  cold  isotherms  corresponding  to 
these  models  are  reviewed  in  the  appendix  and  although  a  similar  study  of  the 
corresponding  temperature  dependent  ones  would  be  of  interest,  its  need  is 
diminished  by  the  recent  appearance  of  a  comprehensive  survev  of  the  thermo¬ 
dynamic  properties  of  matter  at  high  pressures  and  temperatures  by 
g 

S.  G.  Brush.  Extensions  of  the  theory  to  more  complex  systems  is  the  ob¬ 
ject  of  considerable  activity  but  we  cannot  go  here  into  these  recent  develop¬ 
ments  of  the  many -body  theory. 

There  also  exists  many  more  or  less  successful  theories  and  calculations 
of  the  ground  state  energy  and  first  excited  states  of  crystals  near  normal 
densities,  but  little  use  has  been  made  of  these  results  in  semiempirical  EOS 
work,  this  region  being  generally  treated  empirically. 
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On  the  other  hand,  considerable  use  is  made  of  the  semi-classical  statisti¬ 
cal  theories^  (T-F,  T-F-D,  etc. )  in  bridging  the  gap  between  the  experi¬ 
mentally  accessible  region  and  the  simpler  asymptotic  models.  But 

"Although  the  T-F  method  is  known  to  be  approximate,  the  necessary 
analysis  of  the  applicability  of  these  results  has  not  been  carried  out. 

In  the  literature  there  are  only  qualitative  considerations  of  the  non¬ 
applicability  of  the  method  for  small  compressions  (in  the  region  of 
low'  temperatures)  and  on  the  improvement  of  its  applicability  with 
increase  in  temperature.  However,  the  quantitative  problems  of  the 
limits  of  the  regions  of  density  and  temperature  in  which  the  method 
is  applicable  with  a  given  accuracy,  and  on  the  size  of  the  correc¬ 
tions  to  the  quasi-classical  equations  of  state,  remain  essentially 
unresolved. " tt 

This  paper  and  several  others  by  the  same  author  as  well  as  similar  investi¬ 
gations  by  S.  Golden,  N.  L.  Balazs,  Alfred,  etc, ^present  attempts  at  a  solu¬ 
tion  of  this  problem.  In  view  of  the  formal  nature  of  the  developments  used, 
the  conclusions  reached  are  still  tentative  and  the  proper  limits  of  applicability 
of  the  statistical  models  remain  unknown.  The  main  result  reached  in  these 
studies  is  that  inclusion  of  exchange  effects  exactly  within  the  framework  of 
the  semiclassical  model  (T-F-D)  is  inconsistent  and  that  it  would  be  more 
valid  to  keep  exchange  corrections  to  the  lowest  order  since  higher  order 
times  are  of  the  same  order  as  neglected  quantum  corrections.  If  this  is  the 
case,  it  is  questionable  whether  these  statistical  models  have  gr  eater  range  of 
validity  than  their  leading  terms  in  an  asymptotic  expansion  in  inverse  powers  of 
the  volume,  which  can  be  obtained  directly  from  an  electron  gas  model  in  a  suit¬ 
able  background.  It  should  also  be  pointed  out  that  calculations  on  solvable 
models  along  the  statistical  approach  with  quantum  corrections  (harmonic  oscil¬ 
lator,  etc.)  have  been  notably  unsuccessful.  Furthermore  the  conclusions 
reached  in  the  studies  mentioned  above  are  basically  for  isolated  atoms  or  at 
most  for  simple  monatomic  crystals  and  even  less  is  known  about  the  applica¬ 
bility  of  the  statistical  models  to  polyatomic  structures. 

Finally,  many  of  the  semiempirical  EOS  for  solids  are  not  single  phase 
EOS  but  rather  those  of  a  two -component  system:  a  lattice  part  and  an 
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.Extensive  reviews  of  this  approach  are  available  (see  appendix  and  reference  9). 


t:l  IO  1  V  l  t  V  uno  vuvu  u*  v.  muuuiyj 

D.A.  Kirzhnits,  Soviet,  Physic  8,  1081  (1959). 


electronic  component,  the  total  thermodynamic  state  functions  being  obtained 
additively  under  the  assumption  of  complete  thermodynamic  equilibrium. 

Such  a  separation  corresponds  to  the  usual  adiabatic  (or  the  static)  approxi¬ 
mation  in  quantum  mechanics;  in  semiempirical  EOS  work,  the  two  components 
are  generally  taken  to  be  noninteracting,  although  the  Kormer  et  al.  EOS  with 
variable  lattice  specific  heat*1  does  couple  them  empirically.  One  should  note 
that  in  the  high  density  asymptotic  models,  the  separation  is  that  into  a  system 
of  bare  nuclei  and  one  containing  all  the  electrons  while  near  normal  densities 
it  is  usually  one  into  ionic  cores — (nuclei  and  tightly  bound  electrons)  and  a 
system  of  conduction  electrons.  At  the  present  it  is  still  an  open  question 
whether  a  smooth  transition  between  these  two  regions  is  possible  or  whether 
electronic  phase  transitions  (with  or  without  polymorphic  lattice  ones)  must 
necessarily  occur. 

E.  SEMIEMPIRICAL  EOS  FOR  SiOg. 

(1)  In  the  previous  report  (reference  1),  semiempirical  EOS  for  quartz 
and  stishovite  were  given  which  agree  reasonably  well  with  the  experimental 
Hugoniot  data  obtained  in  this  program  and  elsewhere  for  quartz  and  dry 
playas,  if  the  quartz  one  is  used  in  the  low  pressure  region  and  the  stishovite 
one  in  the  high  pressure  region.  Detailed  calculations  for  the  probable  mixed 
phase  region  were  not  made  in  view  of  the  too  great  uncertainties  affecting  the 
quartz -stishovite  transition  (see  section  3) 

The  procedure  adopted  was  the  one  described  in  reference  1: 

(a)  The  cold  isotherm  Ec(V)  was  obtained  in  the  form  of  a  single 
analytical  expression  for  the  whole  range  of  V  ,  the  interpolation  terms  being 
chosen  to  be  compatible  with  the  functional  form  of  the  cold  energy  as  given  by 
the  standard  T-F  asymptotic  model.  The  use  of  an  approximate  analytical 
expression  for  this  asymptotic  (T-F)  Ec(V)  entails  a  negligible  error  as 
compared  to  the  other  uncertainties  of  the  model;  the  largest  one  of  these 
results^  from  the  necessity  of  introducing  an  effective  monatomic  model*  for 
— • 
tsee  reference  1  and  the  appendix. 

♦The  averaging  procedure  adapted  here  differs  from  that  used  by  L.  Knopoff 
and  G.  MacDonald  (Geophys.  Joum.  284,  1958). 
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the  SiO^  complex  if  use  is  to  be  made  of  the  existing  T-F  calculations. 

(Direct  T-F  calculations  for  molecules  are  difficult  and  few  are  available). 

(b)  The  usual  Mie-Griineisen  EOS  with  constant  F  was  adopted  as  a 
thermal  description  of  the  system.  If  no  phase  transition  occurred,  a  value  of 
r  =  2/3  can  be  taken  as  giving  a  crude  behavior  of  SiO^  over  the  whole  V -T 
r°nge  since  T  =  0.65  at  the  reference  state  and  =  2/3  in  the  very  high 
density  and/or  very  high  temperature  asymptotic  models.  A  higher  value  of 
r  =  1.4  for  stishovite  gives  closer  agreement  with  high  pressure  data  for 
quartz  although  any  value  from  2/3  to  1.6  leads  to  calculated  Hugoniots  within 
the  experimental  error  at  these  pressures  ('•500  kbar). 

The  effect  of  a  variable  chemical  composition  (exclusive  of  water  content) 
expressed  as  a  variation  in  effective  Z  and  Vq  is  less  than  that  of  a  change 
in  the  foim  of  the  interpolation  terms  and  of  no  significance  in  view  of  the  over¬ 
all  limitations  of  the  model.  The  low  pressure  cold  isotherm  is  of  course 
sensitive  to  the  bulk  modulus  at  the  reference  state  but  playas  cannot  easily 
be  compacted  to  solid  density  and  no  conclusive  experimental  data  is  avail¬ 
able  on  the  effect  of  small  changes  in  the  chemical  and  mineralogic  composi¬ 
tion  on  the  bulk  modulus  of  the  soild  phase. 

(2)  Other  possible  semiempirical  EOS. 

The  form  of  Ec(V)  used  here  appears  adequate  for  the  present  and  re¬ 
quires  only  three  parameters  determined  from  the  initial  density,  initial 
compressibility  and  sublimation  energy  of  the  solid  (in  addition  to  its  chemical 
composition)  while  converging  to  the  high  density  models*  at  high  pressure. 

A  simpler  form  could  be  derived  in  the  same  way  by  replacing  ihe  E 

c 

T-F  term  by  the  cold  energy  expression  of  an  electron  gas  in  a  positive  lattice 
13  14 

background.  *  This  procedure  would  avoid  the  use  of  an  effective  monatomic 
model  and  allow  a  distinction  between  polymorphs.  The  correlation  pressure 
term  of  this  model  is  known  for  a  few  lattices,  but  the  probably  more  important 
inhomogeneity  correction  is  still  in  doubt.  ^  Further  development  of  such  a 
high  density  model  for  polyatomic  substances  would  be  worthwhile.  The  use  of 


*Nonrelativistic  ones. 
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a  simple  constant  T  Mie-Griineisen  EOS  is  of  course  only  a  crude  first 
approximation.  Several  more  complex  EOS  can  be  used  with  the  same  cold 
isotherm. 

(a)  An  EOS  of  the  form  (4-18),  in  preference  to  (4-15)  for  the  reasons 
given  in  §  III  with,  e.g.  ,t 


G(V,E*)  =  a  + - a— 

cE*p  +1 

where  p  =  V/Vq  and  a,  b,  c  are  adjustable  parameters.  The  parameters 
would  be  determined  as  follows: 

a  +  b  =  T0  (Tat  V  =  Vo,  T  =  0,  P  =  0) 

"a"  chosen  a  priori  to  be  either  2/3  (perfect  gas  asymptotic  limit)  or  deter¬ 
mined  with  c  by  a  least  square  fit  to  experimental  data  and  temperature 
dependent  T-F  calculations  corresponding  to  the  parameters  of  the  effective 
monatomic  Si02  model.  This  empirical  form  would  serve  mainly  to  bring 
the  semiempirical  EOS  thermal  components  closer  to  the  T-F  values  at  high 
temperatures.  If  there  is  no  particular  problem  using  it  for  stishovite  with 
sr  1.4,  it  may  be  troublesome  for  a  -quartz  in  view  of  the  low  initial  T 
value  (0.65)  resulting  from  the  loose  structure  of  this  crystalline  form  of 
Si09  .  This  procedure  would  require  a  small  computer  program  and  in  view  of 

u 

the  relative  insensitivity  of  the  experimental  Hugoniot  to  T  ,  it  should  be 
based  on  a  minimum  of  release  adiabat  points  if  a  greater  range  of  pressures 
and/or  porosity  than  those  examined  are  beyond  the  experimental  possibilities. 

(b)  A  Hildebrandt  EOS  of  the  form  (4-10) 

Assuming  a  purely  temperature  dependent  E*  and  Cv  from  room  tem  - 
peratureup,  one  could  make  use  of  an  EOS  of  the  form  (4-10).  Cy(T)  would 
be  taken  from  experimental  data  at  the  lower  temperature,  and  extrapolated 
with  the  eventual  addition  of  a  purely  temperature  dependent  electronic  com¬ 
ponent.  Such  a  model  could  not  reasonably  be  used  over  a  very  extended 

^Similar  to  Tillotson's  formulation  based  on  (4-15). 
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range  of  compression  if  only  because  of  the  volume  dependence  of  the  electronic 
energy  band  gap  a  id  electron  effective  mass,  but  may  be  of  limited  interest 
for  Hugoniots  up  to  a  few  megabars. 

(c)  Two  components  EOS  of  the  general  form 

dE 

E  =  E  (V)  P  =  -  -r—- 
c  c  c  dV 

eV.T,  -  E*|att+  E*cl  P*(V,  T)  -  P*latt+  P*el 

offer  more  flexibility  in  handling  the  thermal  components  in  semi-empirical 
EOS  work.  The  cold  isotherm  can  be  handled  as  previously  and  separate 
assumptions  can  be  made  on  the  lattice  and  the  electronic  thermal  components 
of  the  state  functions. 

The  lattice  part  can  be  treated  in  any  of  the  ways  described  earlier  (usual 
Mie-Gruneisen,  Hildebrandt  EOS,  etc. ).  In  particular  one  can  use  the  simple 
Dulong-Petit  solid  model  (constant  specific  heat  Cy  =  3R)  from  some  non- 
zero  reference  temperature  Tq  on  'i.e. , 

E  (T)  =  CyT  +  const 

/  E*iatt(To)  \ 

=  3R/T  -  Tq  +  —  3R~  )  =  aR  (T  -  T0> 

P*(V,T)  =  ^3R(T  -  Tq)  =  glatt(V)T  +  const 

For  materials  with  a  relatively  high  Debye  temperature,  such  a  model  is 
rather  unsatisfactory  at  lower  temperature  and  it  would  be  preferable  to  use 
a  temperature  dependent  Cy(T)  ,  dete”mined  from  static  experimental  data 


t 


As  done  by  Altshuler,  et  al  for  metals  ref. 


12,  9. 
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and  a  better  theoretical  model  (Debye  solid,  Debye -Einstein  solid,  etc.).  The 
Hildebrandt  EOS  is  well  adapted  for  not  too  low  temperatures,  giving,  if  Tq 
is  the  reference  state  (room  temperature  e.  g. ) 


P(V,T)  = 


VS'T  (V)  gT  (V) 
o  _  o 
CV(T)  "  Cy(T) 


A 

E‘l»tt<T>  =  /  Cy(T)dt  +  const 


gx  (V) 


"  latt<T>  = 


V 


T  +  const 


The  form  of  Cy(T)  can  be  chosen  to  decrease  at  high  temperatures  from  the 
3R  Dulong-Pctit  value  to  the  3R/2  perfect  liquid  lattice  value*.  Since  such  a 
model  is  not  extendable  to  Tq  =  0  ,  it  is  in  our  opinion  best  to  avoid  the  use 
of  Slater  or  McDugall  type  of  relations  and  treat  STo(v)  empirically. 

The  electronic  thermal  component  can  similarly  be  treated  in  several  ways, 
although  the  quasi-free  electron  gas  model  as  used  by  Altschuler,  et.  al. ,  for 
metals  is  of  course  not  applicable  directly  to  the  semi-conductors  or  insulators 
which  constitute  the  play  as.  Up  to  a  few  ev,  one  can  use  the  various  E*ej , 

P*el  expression.",  developed  for  semi-conductors  e.g. ,  as  (Cornier  et  al,  ** 
or  at  higher  temperatures  the  approximate  analytical  expressions  resulting  from 
a  fit  to  the  T-F  temperature  calculations  of  Latter  and  others,  or  combine  both. 
The  degree  of  validity  of  such  assumptions  is  considerably  more  difficult  to 
assess  than  the  still  unsolved  one  of  the  cold  models,  but  there  is  little  doubt, 
in  view  in  particular  of  the  extensive  Russian  work  on  porous  metals  and  ionic 
crystals,  that  reasonable  EOS  must  include  an  electronic  component  along  the 
above  lines  already  at  medium  pressures  (a  few  megabars),  specially  for 
porous  media. 


*Sec  e.g. ,  Kormer  et  al,  ref.  11. 
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If  a  wide  range  of  porosity  and  pressure  exceeds  the  present  experimental 
possibilities,  it  is  necessary  to  include  in  any  semiempirical  determination  of 
the  thermal  components  of  the  EOS  a  minimum  of  release  adiabat  points,  be¬ 
cause  the  Hugoniot  itself  for  a  solid  or  a  low  porosity  medium  is  relatively  too 
insensitive  to,  e.g. ,  the  I"  chosen  (such  data  were  not  available  in  time  to  be 
included  in  the  theoretical  work  of  this  program). 

F.  CONCLUSIONS  AND  RECOMMENDATIONS 

1.  Almost  any  semiempirical  EOS  formulation  can  be  used  to  represent 
the  Hugoniot  data  in  the  experimentally  accessible  region  and  such  Hugoniot 
points  can  be  reasonably  well  predicted  in  the  absence  of  phase  changes,  from 
a  limited  number  of  parameters  determined  at  the  reference  state. 

2.  Such  low  and  medium  pressure  EOS  can  be  extended  in  a  compatible 
way  to  high  and  very  high  pressures  by  introducing  various  asymptotic  models. 
But,  at  the  present  there  is  no  immediate  way  of  asserting  the  degree  of  validity 
of  the  resulting  EOS  in  any  given  range  of  the  parameters  (outside  the  experi¬ 
mental  region)  and  this  prevents  one  from  selecting  on  sound  grounds  a  particu¬ 
lar  EOS  as  superior  to  others  and  of  making  more  than  qualitative  estimates 
(guesses)  on  their  respective  validity. 

3.  Although  many  so-called  "corrections"  have  been  introduced  with  an 
aim  at  refining  the  original  semiclassical  statistical  T-F  model,  these  are 
often  inconsistent  and  have  not  increased  the  degree  of  reliability  of  the  re¬ 
sulting  modified  statistical  models. 

4.  Recent  developments  in  many-body  quantum  (and  classical)  statistical 
mechanics  give  hope  that  significant  progress  could  be  achieved  now  in  obtain¬ 
ing,  through  bona-fide  estimates,  quantitative  reliability  and  unambiguous  pre¬ 
dictions  at  least  for  the  high  density  models.  Such  theoretical  work,  in  con¬ 
junction  with  a  more  systematic  use  of  existing,  near  normal  density  theories 
and  experimental  work  could  allow  significant  progress  to  be  made  in  resolving 
some  of  the  many  existing  problems  in  EOS  work. 

Before  any  quantitative  reliability  can  be  given  to  proposed  EOS  of  playas, 
it  is  definitively  necessary  to  conduct  a  considerably  more  intensive  investiga¬ 
tion  of  the  eventual  phase  transformations  of  SiO^  and  other  constituent  minerals 
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in  dynamic  compression  and  upon  pressure  release.  In  particular,  any  elec¬ 
tronic  thermal  component  of  the  thermodynamic  state  functions  can  in  no  way 
be  extrapolated  across  such  phase  transitions  from  reference  state  data  on  the 
basis  of  present  knowledge. 

LIST  OF  SYMBOLS 

r  Gruneisen  ratio 
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T  Temperature 

P  Pressure 

E  Internal  energy 
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c 

E*  Thermal  component  of  energy 

p  Density 
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5 -EFFECTS  OF  A  PHASE  TRANSITION  ON  THE  PROPAGATION 
OF  FINITE  AMPLITUDE  WAVES 

G.  E.  Duvall  and  Y.  Horie* 


A.  INTRODUCTION 

Many  materials  have  been  found  to  undergo  a  phase  transition  on  com¬ 
pression.  Under  some  conditions  the  transition  has  important  effects  on  the 
structure  and  propagation  history  of  a  finite  amplitude  stress  wave.  One  such 
effect  is  to  produce  an  instability  in  the  compressive  shock  wave;  another  is 
to  introduce  the  possibility  of  rarefaction  shocks.  In  Section  B  below  the  theory 
of  shook  wave  ‘-lability  is  reviewed  and  extended  to  a  lorm  appropriate  for 
discussion  ol  effects  duo  to  phase  transitions.  In  Sections  C  and  D  polymorphic 
transitions  in  which  the  density  increases  upon  the  increase  of  compression 
are  classified,  their  phase  boundaries  are  described,  and  some  adiabats  and 
R-H  curves  in  the  coexistence  region  are  calculated.  The  theory  is  used  to 
discuss  transitions  in  iron,  bismuth  and  quart/.,  all  of  which  have  been  studied 
experimentally. 

Notation 

P  =  pressure 

V  =  specific  volume 

T  =  temperature  in  degrees  Kelvin 

S  =  specific  entropv 

E  -  specific  internal  energy 

G  =  specific  Gibbs  energy 

H  =  specific  enthalpy 

p  =  mass  density  =  1/V 

F  =  Griineisen  constant 

U  -  shock  propagation  velocity 

u  =  particle  velocity 

*Dept.  of  Physics,  Washington  State  University 
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The  locus  oi  states  in  the  P-V  plane  which  can  be  reached  by  a  single 
shock  from  a  given  initial  state  is  called  the  Hugoniot  or  R-H  curve  for  that 
initial  state. 

The  jump  conditions  for  a  single  shock  running  into  material  in  the  state 
P0  ,  VQ  ,  E0  ,  uQ  =  0  and  compressing  it  to  the  state  P.^  ,  ,  Uj 

are 


U1  =  pl(lJJ  "  ul) 

(5-1) 

1  "  po  P0U1U1 

(5-2) 

=  1<P1 +  po»<vo  -  V1> 

(5-3) 

A  number  in  parentheses  in  the  text  indicates  one  of  the  numbered  references 
at  the  end  of  the  report. 

B.  STABILITY  OF  SHOCK  WAVES 

(1)  Thermodynamic  Criteria 

The  stability  of  a  single  shock  transition  was  first  discussed  by  Rayleigh 
Ref.  (15),  who  concluded  that  only  compressive  shock  waves  were  stable  in 
gases,  "his  conclusion  resulted  from  an  analysis  which  showed  entropy  change 
to  be  positive  for  a  compression  shock  and  negative  for  one  of  rarefaction.  A 
negative  entropy  change  violates  the  second  law,  so  he  concluded  that  only 
compression  shocks  are  stable.  Bethe  Ref.  (16)  investigated  the  stability 
condition  for  a  more  general  equation  of  state  and  found  the  entropy  condition 
to  be  satisfied  if  the  curvature  of  the  adiabat  is  everywhere  positive,  i.  e. 


(5-4) 


He  concluded  that  if  this  condition  were  violated  at  any  point,  that  point  would 
be  one  of  instability  for  break-up  of  the  compression  shock  into  multiple 
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waves.  He  also  demonstrated  that  Condition  (5-4)  can  be  violated  at  a  phase 
boundary,  which  then  emerges  as  a  point  of  instability. 

Minshall  (17)  reported  the  observation  of  multiple  compression  shocks  in 
iron  in  1954,  and  attributed  them  to  plastic  yielding  and  to  a  polymorphic  phase 
transition.  Drummond  (18)  showed  that  violation  of  Condition  (5-4)  leads  to  the 
existence  of  rarefaction  shocks  and  made  calculations  in  iron. 

Rice,  McQueen  and  Walsh  (19)  suggested  a  stability  criterion  which  differ 
from  Cond.  (5-4).  Suppose  OAB  in  Fig.  5-1  is  an  R-H  curve — derived,  per¬ 
haps,  from  experimental  data — and  we  wish  to  test  whether  or  not  a  single 
compression  shock  from  O  to  B  is  stable  against  breakup  into  two  shocks:  one 
from  O  to  A  and  a  second  from  A  to  B.  In  order  to  do  this  we  suppose  the 
single  shock  to  be  unstable,  having  a  compression  profile  like  that  shown  in 
Fig.  5-l(b).  Then  the  velocity  of  the  first  shock  with  respect  to  the  material 
following  is 


(5-5) 


ASSUMED  R-H  CURVE  COMPRESSION  PROFILE  FOR 

ASSUMED  INSTABILITY  AT  A 

FIG.  5-1 
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The  velocity  of  the  second  shock  with  respect  to  the  material  ahead  of  it  is 


(5-6) 


Then  if  -  u^  ,  the  assumption  of  instability  at  A  is  untenable, 

since  the  second  shock  must  overtake  the  first,  in  violation  of  the  assumption. 
Repetition  of  this  test  for  every  point  A  leads  to  the  following  statement:  if, 
for  every  p0<  PA<  PB  . 


(5-7) 


then  a  single  shock  from  O  to  B  is  stable.  Cond.  (5-7)  does  not  appear  to 
be  particularly  useful.  It  is  not  particularly  appropriate  for  a  priori  construc¬ 
tion  of  a  Hugoniot  from  the  equation  of  state;  and  in  normal  practice  it  is  not 
required  for  testing  of  experimental  data,  since  the  conditions  of  the  experi¬ 
ment  reveal  multiple  shock  structures.  If  double  shock  data  were  reduced  by 
single  shock  theory,  Cond.  (5-7)  would  be  useful  in  principle.  In  practice  a 
point  of  possible  instability  would  be  revealed  by  inspection  of  the  P-V  or 

U„  -  u  data. 

S  p 

(2)  Hydrodynamic  Criterion 
In  Reference  (20)  it  is  shown  that  if 


8(u  +  a)/8P 


S 


>  0 


(5-8) 


at  any  point  on  the  path  of  an  accelerating  piston,  then  that  point  may  produce 
a  shock  discontinuity  in  the  flow.  Moreover  this  condition  is  found  to  be  equiv¬ 
alent  to  Cond.  (5-4)  at  each  point  on  the  adiabat  of  the  material. 
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Except  in  very  simple  cases  the  construction  of  a  Hugoniot  requires 
numerical  computations.  This  suggests  the  combined  procedure  for  con¬ 
structing  the  R-H  curve  and  testing  for  instability  to  be  described  in  the  next 
two  sections. 

(3)  Differential  Equation  of  the  R-H  Curve 

We  assume  all  required  equation  of  state  information  to  be  known  and  we 
seek  a  procedure  whereby  the  R-H  curve  can  be  constructed  step-by-step, 
starting  at  the  initial  pressure  and  proceeding  to  higher  values,  testing  at 
each  step  to  determine  whether  or  not  a  single  shock  from  the  initial  state  to 
the  next  higher  pressure  will  be  stable.  The  procedure  for  constructing  the 
R-H  curve  in  this  way  is  available  if  we  have  its  differential  equation  and  a 
suitable  test.  The  differential  equation  can  be  obtained  as  follows: 

Suppose  the  R-H  curve  is  known  up  to  a  pressure  and  that  a  single 
shock  from  (P()  ,  VQ)  to  (P^  ,  V^)  is  stable.  Then  the  internal  energy  at 
<P i  ,  Vj)  is  given  by  the  R-H  equation: 

E!  -  E0  =  !<P1  -  P0><V0  -  Vl>  <5-9' 

If  a  single  shock  to  the  higher  pressure  P^  +  SP^  is  stable,  the  change  in 
internal  energy  will  be,  to  first  order  in  small  quantities, 

SEj  --  jsi'1(v0  -  Vj)  -  i(Pj  +  P^SVj  (5-10) 

A  thermodynamic  path  can  be  found  which  connects  (P^  +  SP^  , 

Vf  +  SV^)  and  (P^  ,  Vj)  .  This  will  require  the  addition  of  heat  and  work. 
The  First  and  Second  laws  of  Thermodynamics  are  to  be  satisfied,  so  that 

dEx  =  TjdSe  -  P1dV1  ,  (5-11) 
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where  any  entropy  change  clue  to  irreversible  internal  processes  is  included 
in  P^dV^  •  Equating  Eqs.  (5-10)  and  (5-11)  yields  a  relation  lor  the  rate  of 
increase  of  entropy  along  the  R-H  curve,  relative  to  that  which  exists  in  the 
adiabat: 


dSe/Dv  =  <V0  -  V1)(dP/dV)R_H/2T1  +  <Pt  -  PQ)/2T1  (5-12) 

if  we  suppose  pressure  to  be  a  function  of  volume  and  entropy,  we  can  express 
the  rate  of  change  of  pressure  with  respect  to  volume  in  any  direction  in  terms 
of  its  partial  derivatives.  In  particular,  for  the  R-H  curve, 


adiab. 


(5-13) 


Eliminating  c!Sc/dV  between  Eqs.  (5-12)  and  (5-13)  yields  the  differential 
equation  of  the  R-H  curve: 


In  obtaining  Eq.  (5-14),  the  following  identities  have  been  used: 

/  8P\  _ / 3P\  /3T\  TT 

\9s/v  \aT/v\as/v  ~  V 


(5-15) 


For  a  material  in  which  compression  is  reversible,  (BP/aV)^  =  (3P/aV)s 

G 

and  S  is  the  total  entropy. 

Integration  of  Eq.  (5-14)  yields  the  R-II  curve,  provided  a  single  shock 
is  stable  everywhere. 
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(4)  Point  Stability  Criterion 


We  wish  now  to  answer  the  following  question:  if  i  single  shock  is  stable 
at  the  pressure  ,  will  a  single  shock  be  stable  to  +  5P^  ?  The 
Rayleigh  line  connecting  (P^  ,  V^)  and  (P^  ,  V^)  intersects  the  R-H  curve 
only  at  these  two  points  (6);  accordingly 


> 

RU,  Px 


(5-16) 


II  P.  is  a  point  of  discontinuity  in  (dP/dV)„  ,  the  value  on  the  lower  side 

1  rv-rx 

of  Pj  is  to  be  taken.  If  the  point  +  8Pj  is  also  to  be  attained  through  a 

single  shock,  Cond.  (5-16)  must  hold  above  ;  i.e. 


-  (—) 
\dV4 


R-H.P 


1 


(5-17) 


Substitution  of  Eq.  (5-11)  into  this  inequality  yields  the  condition  for  stability: 


(Ur  uA)‘ 


l'(V0  -  Vj) 
2V. 


1  "  "  vi> 


>1 


(5-18) 


P  = 


P1  +  SP1 


If 


nv0  -  V1)2V1<  1  , 


(5-19) 


then  Cond.  (5-18)  reduces  to  a  simpler  one:  Pj  ,  V ^  is  a  stable  point  if 


lim 

8l\-H) 

A 


2 

a 


•> 


+  8  Pj 


S  1 


(5-20) 


77 


Since  Cond.  (5-19)  is  normally  satisfied,  Cond.  (5-20)  is  a  useful  one  which 
corresponds  to  the  hydrodynamic  Condition  that  the  flow  be  subsonic  behind  the 
shock. * 

A  situation  in  which  is  a  state  of  instability  is  illustrated  in  Fig.  5-2. 

C.  THERMODYNAMIC  FUNCTIONS  IN  THE  COEXISTENCE  REGION 

We  are  concerned  here  with  instabilities  which  may  arise  when  the 
Hugoniot  curve  intersects  a  phase  boundary  and  the  compressed  material  is 
forced  into  the  coexistence  region  or  on  into  a  second  phase.  Bethe  (16)  has 
shown  that  such  a  point  may  be  a  point  of  shock  instability. 

We  restrict  consideration  to  materials  which  transform  isothermally  to  a 
higher  density  form  when  pressure  is  applied.  Transitions  will  be  classified 
according  to  the  sign  of  dP/dT  in  the  coexistence  or  mixed  phase  region  and 
the  sign  of  dS/dT  on  the  phase  boundary;  procedures  for  constructing 


FIG.  5-2  COMPARISON  OF  SLOPES  AT  A  POINT  OF  INSTABILITY 


♦When  Cond.  (5-19)  is  satisfied,  Cond.  (5-18)  implies  Cond.  (5-20)  and  vice- 
versa.  When  Cond.  (5-19)  is  violated,  a  more  detailed  investigation  based 
on  hydrodynamics  will  be  required  to  establish  the  proper  criteria  for  stability. 
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adiabats  and  R-H  curves  in  the  mixed-phase  region  are  described,  some 
general  properties  are  derived,  and  some  examples  are  treated. 

Although  interest  here  is  centered  primarily  in  polymorphic  transitions  of 
solids,  there  is  no  difference  thermodynamically  between  these  and  the  trans¬ 
formations  of  melting  and  vaporization;  all  are  governed  by  the  Clausius  - 
Clapeyron  equation  in  the  coexistence  region.  This  will  henceforth  be 
designated  "C -region"  or  "C-R",  and  thermodynamic  quantities  in  the  C-region 
will  be  denoted  by  a  subscript  "M".  It  is  also  assumed  that  isothermal  changes 
satisfy  the  Gibbs  condition  at  constant  P  and  T  . 

In  Section  5-1,  a  classification  of  transitions  is  made  according  to  the 
sign  of  dP/ilT  and  dSj/dT  ,  where  the  subscript  "1"  denotes  changes  along 
the  boundary  between  the  less  dense  phase,  denoted  as  phase  one  or  01  ,  and 
the  C-region;  and  relations  between  the  phase  boundary,  the  adiabat  in  01 
and  the  isotherm  in  01  are  derived.  In  Section  5-2  the  slopes  of  adiabat  and 
Hugoniot  in  the  C-region  are  derived  and  some  relations  at  the  boundary  are 
established.  The  slopes  of  adiabats  and  R-H  curves  are  compared  and  the 
curvature  of  the  adiabat  is  discussed. 

For  simplicity  we  assume  the  region  of  interest  to  be  far  from  a  triple¬ 
point. 

(1)  Phase  Boundaries 

The  basis  for  conventional  thermodynamic  treatment  of  phase  transitions 
is  the  assumption  that  extensive  thermodynamic  properties  are  mass-weighted 
averages  of  properties  of  the  two  components,  i.e.,  phases.  Examples  of 
such  are  entropy,  specific  volume,  and  Gibbs  free  energy,  all  of  which 
depend  on  amount  of  material,  and  so  vary  throughout  the  C-region.  Pressure 
and  temperature  can  be  simultaneously  constant  during  a  phase  change  because 
they  are  intensive  variables;  their  values  do  not  depend  on  the  amount  of  ma¬ 
terial  present. 

The  mass-weighted  relation  for  extensive  variables  comes  directly  from 
the  condition  that  P  and  T  are  constant  in  the  C-region  and  that  surface 
energies  can  be  neglected.  Since  dP  =  dT  =  0  ,  then  dG  =  0  =  gjdMj  x 
g  dM.,  ,  where  divl.  is  the  fraction  of  unit  mass  which  goes  into  phase  1  and 

M  M  A 
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dM2  is  the  fraction  which  goes  into  phase  2  in  the  contemplated  change.  But 
-dM^  =  dM2  =  dX  ,  so  g^  =  g2  .  Here  X  denotes  the  fraction  of  mass 
that  has  transformed  from  phase  1  to  phase  2  (Fig.  5-3):  X  =  (V  -  V^)/ 

(V2  ‘  Vl>  * 

With  this  definition  of  X  ,  the  specific  volume  at  point  Q  of  Fig.  5-3  is 


V  =  xv2  +  (1  -  X)Vt 
=  V  +  XAV,  0  Z  X  '<  1 


(5-21) 


where  is  specific  volume  at  point  B  of  Fig.  5-3,  V2  is  its  value  at 
point  E  ,  and 


AV  =  V2  -  V  <  0  (5-22) 

Then  for  constant  T  we  have 

dV  =  AV  dX  ;  dT  =  0  (5-23) 


FIG.  5-3  DEFINITION  OF  THE  TRANSITION 
PARAMETER  X 
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Other  extensive  variables  such  as  H  ,  S  ,  E  ,  are  given  by  similar 
expressions: 


=  Hj  +  X  A 

H 

(5-24) 

=  s  +  x  A 

S 

(5-25) 

-  Ej  t  X  A 

E  . 

(5-26) 

Temperature  on  the  phase  boundary  ABC  of  Fig.  (5-3)  will  be  found  from  the 
Clausius -Clapeyron  equation: 

dP/dT  =  AH/TAV  =  AS/AV  (5-27) 

T  -  Tx  =  f( AV/AS)dP  =  /(dT/dP)dP  (5-28) 

Phase  changes  are  commonly  studied  under  static  or  quasi-static  condi¬ 
tions  in  which  either  P  or  T  is  held  nominally  constant  and  the  other  quantity 
is  slowly  varied  through  the  phase  transition.  The  results  of  such  studies  are 
complete  when  a  curve  P(T)  representing  the  phase  transition  is  obtained, 
along  with  AV(’l')  .  Dynamic  processes  are  usually  other  than  isothermal, 
and  this  fact  leads  to  the  necessity  lor  investigating  the  geometry  of  phase 
boundaries,  coexistence  regions,  adiabats  and  R-H  curves  in  the  S-T  and 
P-V  planes  if  the  relations  between  phase  transitions  and  wave  propagation 
are  to  be  understood. 

We  first  consider  phase  boundaries  in  the  S-T  plane  and  restrict  dis¬ 
cussion  explicitly  to  those  materials  for  which 

AV  <  0  (5-29) 

(dV/9T)p  >  0  (5-30) 
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Here  the  subscript  "P^"  denotes  a  coefficient  at  constant  pressure  in  01 
evaluated  at  the  boundary  between  01  and  the  C-region.  A  derivative  written 
with  a  block  "d"  ,  having  subscript  "1"  on  the  variable  in  the  numerator, 
denotes  a  derivative  along  the  phase  boundary  between  01  and  C-R.  Sub¬ 
scripts  "2"  have  analogous  meanings  for  02  and  its  boundary  with  the  C- 
region. 

The  rate  of  change  of  entropy  along  the  phase  boundary  can  be  expressed 
as 


(5-31) 


(5-32) 


The  transformation  from  Eq.  (5-31)  to  (5-32)  is  made  using  a  Maxwell  rela¬ 
tion  and  the  definition  of  specific  heat  at  constant  pressure.  We  assume  that 


(5-33) 


then,  with  Cond.  (5-30)  we  see  that,  if  dP/dT  <  0  ,  then  dS^/dT  >  0  .  We 
label  this  a  transition  of  Type  1.  If  dP/dT  >  0  ,  then  dS^/dT  may  be 
either  positive  or  negative,  depending  on  the  magnitude  of  dP/dT  .  We  label 
these  transitions  Type  2  and  Type  3,  respectively.  These  classifications  and 
their  properties  are  summarized  in  Table  5-1. 
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Table  5-1 

CLASSIFICATION  OF  POLYMORPHIC  TRANSITIONS 
AV  <  0  ,  (8V/3T)pi  >  0  ,  Cpi  >  0 


TYPE  OF 
TRANSITION 

dP/dT 

AS 

dS:/dT 

1 

<0 

>0 

>0 

2 

>0 

<  0 

>0 

3 

>0 

<0 

<0 

The  slope  ol  the  second  phase  boundary  comes  from  the  identity 


dS  dS  , 

~dT~  =  ~dT  +  * 

if  AH/T  increases  with  T  ,  then  dS2/dT  >  dSj/dT  and  vice-versa.  It  is 
unfortunate  that  knowledge  of  most  phase  transitions  does  not  include  the  sign 
and  magnitude  of  this  derivative.  It  will  be  neglected  in  most  of  what  follows. 

From  Table  5-1  it  is  possible  to  construct  schematic  representations  of 
the  geometry  of  the  phase  boundaries  in  the  S-T  plane.  These  are  shown  in 
Fig.  5-4.  For  Type  1,  Fig.  5 -4(a)  0 2  lies  above  and  to  the  left  of  01  since 
AS  >  0  ;  consequently  temperature  decreases  along  an  adiabat  AC  traversed 
from  01  to  02  .  Type  2,  shown  in  Fig.  5 -4(b)  has  the  relative  locations  of 
01  and  02  interchanged,  AS  <  0  and  T  increases  from  01  to  02  along 
an  adiabat.  Type  3,  Fig.  5 -4(c)  has  phase  boundaries  of  negative  slope,  01 
is  above  and  to  the  right  of  02  ,  and  T  decreases  from  01  to  02  along  an 
adiabat. 

Bridgman  (21)  has  noted  that  dP/dT  for  polymorphic  transitions  is 
commonly  the  order  of  50  oars  per  degree  Centigrade;  i.e.  transistions  of 
Types  2  and  3  are  more  common  than  those  of  Type  1.  Slater  (22)  points  out 
that  this  is  to  be  expected  since  increases  in  entropy  are  normally  associated 
with  increases  in  volume.  He  goes  farther  to  show  that,  if  AS/AV  is 
identified  with  (8S/8V)T  for  a  single  phase,  we  have 

(8S/8V)t  =  TCV/V  (5-34) 


83 


where  F  is  the  Griineisen  constant  and  is  specific  heat  at  constant  volume. 
For  Cy  -  3R  per  mol  ,  T  =  2  and  V  =  10  cc/mol  , 

AS/AV  =  (8S/8V)t  ^  50  bars/°C  (5-35) 


Using  Eq.  (5-34)  in  Eq.  (5-32)  and  equating  Cp  to  ,  we  arrive  at  the 
expression 


cp  tt  /av\ 

t  [  v  larL 


(5-36) 


Taking  F  =  2  and  (l/V)(8V/8T)p  *<  3  x  lO-^  ,  we  conclude  that 

dSj/dT  ^(C  /T)|l  -  0.6  x  10'3  T  ]  (5-37) 


for  most  solids.  Thus,  except  for  very  high  temperature  transitions,  we  will 
normally  expect  Type  2  to  prevail.  In  fact  no  transition  of  Type  3  is  known  to 
the  authors. 

Incidentally,  the  above  classification  suggests  a  fourth  Type  in  which 
dS/dT  is  negative  with  02  above  01  ;  this  corresponds  to  AV  >  0  ,  which 
is  a  case  of  no  interest  at  present. 

A  procedure  for  determining  the  relative  dispositions  of  the  pure  phase 
regions  and  their  boundaries  in  the  P-V  plane  is  suggested  by  a  calculation 
in  Reference  (16).  Considering  specific  volume  to  be  a  function  of  P  and  T, 
we  nave 


(5-38) 


85 


+ 


(5-39) 


av\  dT 
9T/pidP 


For  (3T/3P)S1  , 


(5-40) 


The  inequality  follows  from  Conds.  (5-30)  and  (5-33)  and  implies  that 

(8V/3P)S1  >  OV/3P)T1  (5-41) 


Conversely,  this  inequality  implies  that  of  Eq.  (5-40).  Eqs,  (5-38),  (5-39) 
and  (5-40)  then  imply  the  following  relations  among  dV^/dP  ,  (9V/9P)S1  , 
and  (8V/9P)T1  for  the  three  types  of  phase  transitions  described  in  Table  5-1; 

Type  1 


dP 

dT 


<  0 


(5-42) 


Type  2 


.  ,  dT/dP  >  (3T/8P)gl  and  (5-43) 
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Type  3 


3v\  dP 
3TypidT 

(5-44) 


These  relations  are  described  in  Fig.  5-5.  For  Type  2,  dV^/dP  may  be 
positive;  for  the  others  it  must  be  negative  if  (3V/3P)S  and  (3V/3P)^,  are 
both  negative. 

The  slope  of  the  second  phase  boundary  is  obtained  from  the  definition 

V„  -  V  =  AV  =  f(P)  <  0  .  Then 
^  1 

dV2/dP  -  dV^/dP  +  d(AV)/dP 


1  + 


dPj/dV 

dPl  d(AV) 
dV  dP 


(5-45) 


if  d(AV)/dP  >  0  and  (dPj/dV)  |d(AV)/dp]  >  -  1  ,  then  dP2/dV  <  0  and 
IdP^/dVl  >  |dPj/dV|  .  This  follows  from  the  previous  result  that 
dPj/dV  <  0  when  dP/dT  <  0  .  However,  if  (dPj/dV)  |d(AV)/dpj  <  -  1  , 

then  the  denominator  of  Eq.  (5-45)  is  negative  and  dP2/dV  >  0  .  This  be¬ 
havior  would  be  expected  near  a  critical  point  where  AV — >0  as  P — ^P^  . 

For  purposes  of  illustration  we  consider  a  region  where  dP2/dV  = 
dPj/dV  .  Then,  for  example,  for  Type  1  there  may  be  a  situation  like  that 
shown  in  Fig.  5-6.  Here  part  of  the  region  contained  between  the  two  phase 
boundaries  is  triplymapped:  first  by  the  isotherm  GB  of  phase  1,  then  by 
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TYPE  I 


0 

OV/3P)s| 

ov/aP)TI 

dV,/dP 


0 

dV(/dP 

(av/dP)51 


TYPE  2 


OV/3P)T| 


TYPE  3 


o 

(dV/9P)s, 

dV,/dP 

OV/SP)T| 


FIG.  5-5  RELATIVE  SCOPES  OF  ISOTHERM,  PHASE  BOUNDARY 
AND  ADIABAT  IN  P-V  PLANE 
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the  isotherms  BE  of  the 
mixed  phase,  and  lastly  by  the 
isotherms  EH  of  Phase  2. 

EH  must  lie  above  DEF  if 
OV/8T)p2  >  0  and  Cp2  >  0 
by  virtue  of  the  equation 
corresponding  to  Eq.  (5-39) 
for  the  second  phase  boun¬ 
dary  when  dP/dT  <  0  . 

Any  apparent  anomaly  in 
such  a  multiple  mapping  of  the 
P-V  plane  is  immediately 
resolved  by  recalling  that 
behavior  in  the  P-V  plane  simply  reflects  a  projection  of  relations  on  the 
P-V-T  surface.  Fig.  5-7  represents  a  P-V-T  surface  for  a  material  with 
a  phase  change  of  Type  1.  AB  is  the  phase  line  in  the  P-T  plane, 
dP/dT  <  0  ;  it  is  the  projection  of  the  cylindrical  surface  of  the  coexistence 
region,  DFB'B",  onto  the  P-T  plane,  B"D'  is  the  projection  of  B"D  onto  the 
P-V  plane,  etc. 

The  three  types  are  summarized  graphically  in  the  P-V  and  S-T  planes 
in  Fig.  5-8. 


FIG.  5-6  TRIPLY-MAPPED  REGION  IN  THE 
P-V  PLANE  FOR  TYPE  1 
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FIG.  5-7  EQUATION  OF  STATE  SURFACE  IN  P-V-T  SPACE 
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ISOTHERM 


MULTIPLY 
MAPPED  REGION 


REGION 


ISOTHERM 


ADIABAT 


ISOTHERM 


(2)  Adiabats  and  R-H  Curve 

(2.  1)  Adiabat.  We  first  obtain  :»n 
expression  for  the  adiabat  traversing  the 
coexistence  region  in  the  P-V  plane.  Re¬ 
ferring  to  Fig.  5-9,  we  wish  to  calculate  the 
volume  at  point  C,  pressure  P,  which  has 
the  same  entropy  as  point  A  on  the  phase 
boundar> .  In  order  to  do  this  we  compute 
the  entropy  at  C  by  integrating  along  ABC 
and  setting  S  -  SA  =  0  .  The  resulting 


FIG.  5-9  ADIABAT  IN  THE 

COEXISTENCE  REGION 


equation  can  be  differentiated  and  solved  for  the  adiabatic  slope  at  C  ;  i.e. 


0  =J  (dS1/dP)dP  +  XAS 


(5-46) 


where  AS  =  S0  -  S,  evaluated  at  P  =  PD  .  Henceforth  the  quantities  at 

Si  D 

B  and  B'  are  denoted  by  subscripts  "i"  and  "2";  quantities  at  the  bountlary 
point  A  will  be  denoted  by  subscript  "A". 

From  Eq.  (5-21) 


X  =  X(P)  =  (V  -  y})/(y2  -  Vx)  (5-47) 


Combining  Eqs.  (5-46)  and  (5-47): 


B 

(dSj/dPJdP  +  (AS/AV)(V  -  V  )  =  0  (5-48) 


Differentiate  this  to  obtain 


tlSl  AS  /8V\ 
dP  Xvlap/SM 


dV 

dP 


+  (v  _  v 

'v  vrdP\AV/ 


o . 


(5-49) 


92 


Here,  as  before,  total  derivatives  refer  to  variations  along  the  phase  boundary 
adjacent  to  phase  1.  All  quantities  are  functions  of  the  pressure  P,  and  the 
subscript  SM  means  entropy  is  constant  and  the  quantity  is  evaluated  in  the 
mixed  phase  region.  Recall  that  AS/AV  =  dP/dT  and  solve  for  the  adiabatic 
derivative: 


dVl  ^fldT  v  . _d_  /.  dP\ 

dP  dP  dP  ‘  (V  "  Vl)dP  (  n  dT/ 


(5-50) 


This  is  the  required  adiabatic  slope.  The  derivatives  dV^/dP  and  dS^/dP 
are  given  by  Eqs.  (5-39)  and  (5-51)  following: 


dSl  _  /_3S  \  /_3S_\  dT 

dP  ’  \  3P/T1  \  9T/pi  dP 


(5-51)  ^ 


/  3V\  Si  dT 
\  3T/pi  T  dP  * 


(5-52) 


Substitution  of  Eqs.  (5-39)  and  (5-52)  into  Eq.  (5-50)  yields 


/  3V\ 

/  3V\ 

q/3V  ' 

\  dT  CP1 

/dT\2 

T  If  \  d 

UpL. 

\3pL, 

2\3T, 

nidP  T 

(dP/ 

-  <v  -  vl>dP 

(ln  dT/ 

(5-53) 


Eq.  (5-53)  applies  to  an  arbitrary  point  in  the  coexistence  region.  At  the 
boundary  of  phase  1,  V  =  ,  and  it  is  seen  that  a  discontinuity  in  the 

adiabatic  slope  exists,  as  shown  by  Bethe  (16);  i.e,  at  the  phase  boundary 


(5-54) 
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J 


% 


(5-55) 

where  the  subscript  SMI  refers  to  the  mixed  phase  side  at  the  phase  boundary. 
Subtracting  the  first  of  these  from  the  second  yields 


T  /  9V\  2  /aV\  dT  CPl/dT\  2 

cpil  8TyPl  '  aTyP!dP  T  'dP' 


(5-56) 


PI  dT 
T  dP 


2 

>  0  . 


Consequently  there  is  always  a  di  continuity  in  the  adiabatic  slope  at  the  phase 
boundary,  except  when  (dT/dP)gl  =  dT/dP  ;  and  the  slope  of  the  adia'oat  in 
the  pure  phase  is  always  steeper  than  that  in  the  mixed  phase,  i.  e. 


9P/  8v|  gl  > 


|ap/avl 


SMI 


(5-57) 


The  relation  of  (9V/9P)g^  to  other  slopes  at  the  phase  boundary  can  be 
determined  for  the  three  Types  of  transition  as  follows: 


Type  1 

According  to  Eqs.  (5-39)  and  (5-54), 


/_2Y|  =  /3V\  il  £pi/JT\  2 

aP'sM!  "  dP  '8T/pidP  '  T  \dP/ 


(5-58) 
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Type  2 


Eq.  (5-54)  can  be  written 


(5-59) 


By  the  argument  leading  to  Eq.  (5-37),  we  may  expect  that  the  right  hand  side 
of  Eq.  (5-59)  is  normally  negative,  though  it  may  have  either  sign. 

Type  3 


(5-60) 


The  relations  between  various  slopes  at  the  phase  boundary  are  displayed  in 
Fig.  5-10. 

The  curvature  of  the  adiabat  can  be  obtained  by  differentiating  Eq.  (5-50) 
directly: 


_d 

dP 


dV. 


dP 


dSldT 
dP  dP 


(V  -  v2) 


(5-60.1) 


At  the  phase  boundary  where  V  =  ,  this  reduces  to 


_d 

dP 


dV. 


dP 


dSl  dT 
dP  dP 


(5-60.2) 
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0 


TYPE  I 


(3V/9P)SI 

OV/9P)T| 

dV,  /dP 


(9V/9P) 


SMI 


TYPE  3 


o 

(9V/9P)j, 

OV/3P)SM1 

dV,/dP 

(9V/9P)ti 


FIG.  5-10  RELATIVE  SLOPES  AT  THE  PHASE  BOUNDARY 
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(2.2)  R-H  Curve.  The  R-H  curve  and  adiabat  in  a  normal  materia] 
of  single  phase  have  a  second-order  contact  at  the  origin  of  the  Hugoniot  (9); 
the  entropy  difference  between  Hugoniot  and  adiabat  is  then  third  order  in  the 
compression.  If  the  intersection  of  phase  boundary  and  Hugoniot  is  a  point  of 
instability  for  a  single  shock  with  final  amplitude  in  the  C-region  or  phase  2, 
then  this  intersection,  A  in  Fig.  5-11,  may  be  expected  to  serve  as  initial 
state  for  a  second  shock  following  the  first  (j/g  in  Fig.  5-11).  Accordingly  the 
relation  between  that  portion  of  the  Hugoniot  originating  at  point  A  and  the 
adiabat  in  the  C-region  originating  at  A  may  be  expected  to  be  the  same  as 
that  existing  between  adiabat  and  Hugoniot  in  a  single  phase.  A  careful  review 
of  the  premises  on  which  the  single -phase  result  is  based  reveals  no  reason 
for  doubting  that  it  holds  for  the  C-region,  and  a  direct  calculation,  based  on 
the  following  analysis,  verifies  it. 

In  order  to  determine  the  locus  of  the  R-H  curve  in  the  C-region,  consider 
the  situation  depicted  in  Fig.  5-12.  Point  B  represents  the  final  state  of  a 
second  shock  originating  at  point  A  on  the  phase  boundary.  The  enthalpies  at 
A  and  B  are  connected  by  the  R-H  equation;  they  can  also  be  calculated  by  an 
integration  along  the  path  ADB.  Equating  the  two  yields  the  equation  of  the 


v 


GA -  5039-32 


GA-S098-SS 


FIG.  5-11  SECOND  SHOCK  ORIGINATING 
AT  A  PHASE  BOUNDARY 


FIG.  5-12  CALCULATION  OF 
HUGONIOT  IN 
C-REGION 
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R-H  curve  in  the  C -region: 


F 

HB  -  ha  =  “2“ (p  -  PA)(V  +VA)  =  f  (dHj/dPKIP  +  XAH  .  (5-61) 

•'PA 


where  H^(P)  is  enthalpy  on  the  phase  boundary  and 


V  -  V 

v  v  1  HP 

AH  ■  T<v  -  vi»-3?- 


(5-62) 


The  right  hand  side  of  Eq.  (5-62)  is  a  function  of  P  alone.  Substituting 
Eq.  (5-62)  into  Eq.  (5-61)  and  differentiating  with  respect  to  P  yields  an 
equation  for  the  slope  of  the  Hugoniot  in  the  C -region. 


/dv\ 
idPi 


RH 


dP  1 

m  -  w  -  pA> 


m  dV.  dH 
T  dP  _ 1  _ 1 

dT  dP  "  dP 


T<V-Vl»fp(i)+2<VA-V>  +  Vl 


Divide  Eq.  (5-63)  by  T  dP/dT  and  substitute  into  it  the  identity 


(5-63) 


dHj/dP  =  V  +  T  dS1/dP 


(5-64) 


to  obtain 


¥) 


RH 


2T  dP  (P  '  PA)\dP/RH 


^dv\ 


dV,  dS 


dP 


__1  dT 
dP  dP 


,  d  /,„  dP\  ,  VA  "  V  dT 
(V  '  Vl)  dP  (  n  dT  /  2T  3P 


(5-65) 
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The  slope  at  the  phase  boundary,  (dV/dP)DU1  ,  is  given  by  Eq.  (5-54)  or  by 

It  ii  1 

setting  V  =  in  Eq.  (5-65);  the  curvature  at  the  phase  boundary  is 

given  by  Eq.  (5-60.2). 

(2.3)  Relative  Slopes  of  Adiabats  and  R-H  Curve.  For  material  in 
a  single  phase,  the  slope  of  the  Hugoniot  and  the  slope  of  an  adiabat  crossing 
the  Hugoniot  are  related  by  Eq.  (5-14),  which  can  be  cast  in  the  following  form 

1  ‘  *U  ^  ’  (5-66) 

c 


■  2V<V0-V> 


Here 


c2  =  -V2(dP/dV)RH 
a2  =  -V2(dP/8V)s 


r  =  (v/cv)op/3T)v  . 


A  relation  identical  in  form  can  be  derived  for  the  C -region.  Eq.  (5-65)  can 
be  written: 


1  _  —  —  /p  -  p  ) 

2T  dP  '  A' 


1  + 


(5-67) 


where  (P^  ,  V^)  represents  the  intersection  of  the  01  phase  boundary  and 
the  Hugoniot. 
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2  2 
Define  c  and  a 


as  in  Eq.  (5-66),  then  Eq.  (5-67)  becomes 


2 

c 


rM 


2V  <VA  -  V> 


(5-68) 


where 


v  dp  _  v/dp\  dT  T/as\ 

CVMdT  Tt»V/sMdP  VM  WvM 


/a p\  .  /Vp\  .  a p 

13T/VM  \  /  SM  '  dT  • 


(5-69) 


The  specific  heat  in  the  C-region  can  also  be  expressed  in  the  form: 


dsi  d2p 

CVM  ■  TdT  +  T<V  -  Vl>^3 

dl 


-  T 


(5-70) 


This  is  obtained  by  differentiating  Eq.  (5-46). 

2  2 

In  order  to  determine  the  sign  of  c  -a  in  Eq.  (5-68),  note  that  the  sign 

of  the  bracket  is  positive  if  the  slope  of  the  R-H  curve  at  P  is  greater  than 

2  2 

that  of  ihe  Rayleigh  line  from  .  This  is  certainly  true  if  (d  P/dV  )j^j^m  >0; 

this  condition  is  unnecessarily  restrictive,  but  we  shall  use  it  for  simplicity. 

2  2 

Similarly  if  (d  P/dV  <  0  everywhere  between  PA  and  P  ,  the 

bracket  is  negative.  In  order  to  determine  the  sign  of  -(dV/dP)g^  ,  refer 
to  Eq.  (5-50);  it  can  be  written 


-  <V  -  VI>dP 


(5-71) 
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Here  0V/3P)g^j ^  is  evaluated  at  V^(P)  ;  it  is  not  the  slope  of  the  adiabat  of 

Eq.  (5-71)  where  it  intersects  the  phase  boundary.  Then  if  (8V/8P)C  <0 

bl 

at  all  points  of  the  phase  boundary,  so  is  (BV/SPJg^n  according  to  Fig.  5-10. 
Then  if  the  second  term  on  the  right-hand  side  of  Eq.  (5-71)  does  not  override 
the  first,  (8V/8P)SM  <  0;  VA  ”  V  is  Pos^ive*  80  that 

sgn  (c2-a2)  =  sgn  (dT/dP)  sgn  (d2P/dV2)RHM  »  (5-72) 


where 


sgn  x  =  -t-  lifx  >0 
sgn  x=-lifx<0. 


For  «l2P/dV2>RHM 


>0  : 


Type  1 , 


c 


2 


<0 


Types  2  and  3, 


2  2 

The  relation  between  c  and  a  is  normal  for  Types  2  and  3,  i.e.  the  same 

2  2 

as  for  a  single  phase  with  (0  P/8V  )g  >  0  .  Type  1  is  abnormal,  but 
corresponds  to  the  inversion  of  the  order  of  adiabats  in  the  C -Region.  This 
is  illustrated  in  Fig.  5-13.  A  single  shock  from  point  A  into  the  C-region  is 
stable.  Type  2  exhibits  the  same  structure  as  a  normal  single  phase  material. 
Type  3  is  interesting  because  the  phase  boundary  splits  the  isotherms  and 
adiabats  (Cf.  Fig.  5-10).  However,  the  geometric  structure  is  not  violated 
by  the  assumption  of  positive  curvature  and  Cond.  (5-72).  This  is  illustrated 
in  Fig.  5-14. 
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FIG.  5-13  ANOMALOUS  R-H  CURVE  FOR  TYPE  1  TRANSITION 

2  2  2  2 
The  assumption  that  (d  P/dV  )RHM  <  0  reverses  the  sign  of  c  -  a 

in  Cond.  (5-72);  this  leads  to  instability  in  the  shock  for  all  Types. 

This  calculation  can  be  verified  in  the  neighborhood  of  the  boundary  by 
expanding  V(S,  P)  in  a  series  about  the  boundary  point  A: 
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AOIABATS 


FIG.  5-14  ADIABATS  AND  HUGONiOT  FOR  TYPE  3 
TRANSITION 


S  -  can  be  expanded  in  powers  of  (P  -  P^)  along  the  Hugoniot  in  the 
manner  described  by  Courant  -  Friedrichs  (23)  with  the  result: 


S 


— 1 —  iLX  /p  _  p  )**  + 
12T .  _2  {y 


A  dP 


(5-74) 


Entropy  increases  along  this  segment  of  the  Hugoniot  as  P  increases,  for 
d2V/dP2  >  0  . 

From  Maxwell's  relations, 

<dV/8S)p  =  OT/9P)s  =  dT/dP  (5-75) 
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in  the  C -region.  Consequently  Eq.  (5-73)  becomes 


1 

12T. 

A 


dT  d2V  \ 

dP-VA 


(p 


(5-76) 


Differentiating  Eq.  (5-76)  yields 


RHM 


(5-77) 


2  2 

For  (d  V/dP  k  >  0  ,  the  right  hand  side  of  Eq.  (5-77)  has  the  sign  of 
dT/dP  .  This  is  in  accord  with  Cond.  (5-72)  since  (dV/dP)^™-  - 
OV/dP)SM  has  the  same  sign  as  c  -  a  . 


(2. 4)  Curvature  of  Adiabats.  The  curvature  of  the  adiabat  in  the 
C-region  depends  upon  second  order  thermodynamic  coefficients,  and  it  is, 
accordingly,  difficult  to  make  any  general  statements  about  its  sign.  Eq.  (5-56) 
can  be  written 


(5-78) 


so  we  can  say  that  if  (3V/9P)C.  increases  with  increasing  pressure  along  the 

^  2 

phase  boundary,  and  if  T(dS^/dP)  /Cpj  does  not  increase  more  rapidly  than 
(dV/8P)sl  ,  then  (3V/3P)SM1  increases  with  P  so  (82V/aP2)SM1  >0* 

Eq.  (5-60. 1)  provides  an  explicit  statement  of  the  curvature.  It  is  of  little 
help,  however,  for  the  same  reason  that  affects  Eq.  '5-78):  d  S^/dP  is  not 
normally  known.  An  expression  for  the  curvature  in  terms  of  more  directly 
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observable  quantities  can  be  written  as  follows: 


+ 


-  2 


"PI 


dT  d2T 
dPdP2 


(5-79) 


+ 


A 

dP 


(V  -  vx) 


This  equation,  too,  contains  too  many  undetermined  coefficients,  although 
they  are  in  principle  known.  Bethe  (16)  has  examined  the  separate  terms  of 
Eq.  (5-79)  for  some  known  phase  transitions,  and  he  speculates  that  the 
curvature  is  always  positive.  However,  no  such  firm  conclusion  is  justified 
without  more  exhaustive  study. 


D.  EXAMPLES 


Some  data  pertaining  to  phase  transitions  in  bismuth  (BiMBill)  ,  iron 

(a— h.  c.p. )  and  quartz  (a  —  stishovite)  are  given  in  Table  5 -II.  Subscript 

"A"  refers  to  the  intersection  of  the  phase  boundary  and  the  Hugoniot  with 

origin  at  room  temperature  and  atmospheric  pressure.  Values  of  specific 

heat  and  thermal  expansion  coefficient  are  taken  at  atmospheric  pressure. 

Values  of  dVj/dP  are  determined  from  Eq.  (5-39);  values  of  (dP/dT)^^ 

are  obtained  from  measured  values  of  the  R-H  curve  above  the  phase  boundary 

24 

in  the  manner  described  by  Duff  and  Minshall.  The  procedure  is  to  extrapo¬ 
late  the  R-H  data  to  the  phase  boundary  and  from  this  determine  (3V/8P)gMi  s 
(dV/dP)DU1  .  Substitution  of  this  into  Eq.  (5-54)  along  with  the  other  thermo- 
dynamic  data  makes  it  possible  to  calculate  (dP/dT)  s  (dP/dT)cajc  .  A 
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useful  form  of  Eq.  (5-54)  is  that  given  by  Duff  and  Minshall: 


/dP\2  2q1  dP  CP1 

ldT/  *RH1  ‘  *1  dT  '  TV<0RH1  -  ^ 


=  O 


(5-80) 


where  the  following  abbreviations  have  been  used: 


o  =  _1  /dV\ 

^RHl  -  -  V^dP^ 


/3  =  -±  tm 

V1  \3P/t1 


(5-81) 


=  _L  (SOL) 

V  '3T/pi  ’ 


By  substituting  directly  measured  values  of  dP/dT  into  ]  (5-54)  it  is 

possible  to  calculate  0ADM1  =  (l/V)(9V/3P)gM1  .  represents 

equilibrium,  then  it  should  be  true  that  =  ^DMl  *  The  failure  of 

this  inequality  to  hold  reflects  the  difference  between  dP/dT  and  (dP/dT)  j . 
These  differences  may  result  from  nonequilibrium  effects  in  shock  compression, 
from  errors  in  extrapolation  of  the  shock  data  back  to  the  phase  boundary,  or 
from  nonhydrostatic  stress  distribution  in  shocked  material.  Construction  of 
the  equilibrium  R-H  curve  in  the  C -region  can  be  accomplished  by  substituting 
the  appropriate  thermodynamic  coefficients  into  Eq.  (5-61);  an  alternative 
procedure  is  to  examine  Eq.  (5-65)  for  deviations  from  constant  slope. 
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TABLE  5 -II 

PHYSICAL  DATA  ON  PHASE  TRANSITIONS 


Bismuth 

Iron 

Otz-Stishovite 

V  ,  cc/g 

.  1020 

.  128 

.3775 

Cpi,  kb  cc/g°C 

1.  21  x  10"3  (10) 

4. 15  x  10-3 

9.2  x  10"3 

-6 

-6 

-6 

°C 

40.  2  x  10  (10) 

36.3  x  10 

38.4  x  10 

f\,  kb'1 

2.46  x  10"3  (10) 

0.51  x  10"3  (11) 

0.76  x  10"3 

PA*  kb 

25. 3* 

132  (11) 

144  (14) 

VA,  cc/g 

.0929  (13) 

.1197  (11) 

.3145  (14) 

TA,  *k 

315  (10) 

310  (11) 

476  (14) 

dP/dT,  kb/°C 

-.0500,  -.0508 

-.  075  (12) 

+.0177  (15) 

AV,  cc/g 

-.  0047 

-.0041  (11) 

(-1/V)  dVj/dP 

2.86  x  10"3 

.99  x  10"3 

-1.5  x  10‘3  : 

0RH1’  kb 

13  x  10  "3 

2. 18  x  10~3 

1.59  x  10"3 

^ADMl’  kb 

21  x  10"3 

23  x  10 '3 

194  x  10'3 

<dP/dT>calc'  kb/°C 

-.067 

-.29 

+.225 

♦Adjusted  to  static  value. 

Numbers  in  parentheses  indicate  references  at  end  of  paper. 

If  no  reference  is  indicated,  numbers  are  from  Smithsonian  tables  or 
calculated  here. 

=  isobar ic  volume  expansion  coefficient. 

=  isothermal  compressibility  of  pure  phase. 

=  adiabatic  compressibility  of  mixed  phase. 

=  (-l/V)(dV/dP)RH1 


*1 

^ADMl 

^RHl 
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(1)  Bismuth 


24 

Measurements  reported  by  Duff  and  Minshall  were  assumed  to  represent 
the  Bil-Bill  transition  and  were  of  two  kinds.  They  varied  initial  temperature 
of  the  bismuth  specimen  and  measured  the  pressure  amplitude  of  the  first 
wave.  From  their  data  they  inferred  a  value  of  -0.  0508  kbar/°C  for  the 
dP/dT.  This  compares  very  favorably  with  Bridgman's  value  of  -0.0500.  The 
amplitude  of  the  first  wave  was  about  3.  5  kbar  greater  than  the  pressure  re¬ 
ported  by  Bridgman  for  the  same  temperature.  This  difference  might  be  due  to 
a  rate  effect  in  the  transition  or  to  nonhydrostatic  stress  or  to  both.  A  limited 
attempt  revealed  no  decay  of  the  transition  pressure  with  travel  distance, 
which  would  indicate  the  absence  of  any  but  slow  rate  effects. 


Second,  they  measured  the  amplitude  of  the  second  shock  for  different 
shots  with  the  same  initial  temperature  and  extrapolated  the  results  to  obtain 
»  given  in  Table  5 -II.  Extrapolation  errors  are  estimated  at  60  percent 

Knl 

by  Duff  and  Minshall. 


With  the  data  of  Table  5-II  the  change  of  (dV/dP)Rp  within  the  C -region 
can  be  estimated.  From  the  given  value  of  $ADM1  we  find  at  the  boundary 
of  phase  1: 


RH1 


1.95  x  10 


cc/g  kb 


Assume  that  Cpi  =  CV1  ,  dP/dT  =  const.,  dV^/dP  =  const,  and  that 
,  Cp^  are  constant.  Since  dP/dT  =  const., 

(T  -  Ta)  =  -  20  (P  -  PA)  . 

Integration  of  Eq.  (5-65)  with  the  above  assumptions  yields  the  solid  curve  of 
Fig.  5-15.  The  difference  between  this  and  the  measured  R-H  curve  is  notable 
because  of  its  magnitude  and  because  the  curvature  of  the  computed  curve  is 
negative.  The  latter  result  may  well  be  due  to  neglect  of  variations  of  thermo¬ 
dynamic  coefficients.  The  former  is  so  large  that  one  is  almost  forced  to  con¬ 
clude  that  kinetic  effects  are  entering  into  the  shock  compression,  particularly 
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SPECIFIC  VOLUME, V— cm3/g 


FIG.  5-15  COEXISTENCE  REGION  FOR  Bil-Bill 

when  it  is  recalled  that  the  pressure  of  the  phase  boundary  measured  in  shock 
compression  exceeds  that  reported  by  Bridgman. 

A  comparison  of  static  and  dynamic  data  near  the  transition  pressure  is 
shown  on  an  expanded  scale  in  Fig.  5-16. 

(2)  Iron 

Minshall  described  a  double  shock  structure  in  iron  in  1954  (17)  and  suggested 
it  might  be  due  to  a  phase  transition  at  130  kb.  These  and  other  measurements 
were  later  reported  in  detail  by  Bancroft  and  others  (25)  and  analyzed  in  more  de¬ 
tail.  Johnson,  Stein,  and  Davis  reported  an  extensive  setof  measurements  in  1961 
(26)  in  which  initial  temperature  was  varied  from  liquid  nitrogen  to  nearly  the 
a  -  y  temperature.  They  suggested  on  this  basis  that  the  130 -kb  transition  was  a 
transition  to  a  new  phase,  and  this  has  since  been  verified  by  x-ray  diffraction 
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SPECIFIC  VOLUME, V - cm5/j 

• A- SOI#- 87 

FIG.  5-16  PHASE  DIAGRAM  FOR  BISMUTH 

measurements  at  static  high  pressure  and  the  new  phase  is  found  to  be  h.  c.  p. 
From  their  data,  dP/dT  =  -0.  075  kb/°C.  Substitution  of  this  value  into 
Eq.  (5-80)  gives  0ADM1  shown  in  Table  5-n. 

Minshall's  data  include  measurements  of  the  second  shock  amplitude. 
These  can  be  extrapolated  to  obtain  as  reported  in  Table  5-11.  This 

is  only  10  percent  of  0ADM1  ;  the  difference  is  much  greater  than  was 
found  for  Bi  .  The  value  of  (dP/dT)calc  derived  from  flRH1  differs 
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correspondingly  from  the  directly  measured  value.  Remarkably  enough  it  is 
the  same  as  dP/dT  for  the  a  -  y  transition  as  determined  by  Claussen  and 
others  (26).  The  phase  diagram  for  iron  is  shown  in  Fig.  5-17. 

It  is  difficult  to  attribute  such  differences  to  rate  effects  in  shock,  since 
the  first -wave  amplitudes  agree  very  well  with  the  static  values.  Moreover 
there  is  little  evidence  for  decay  of  the  130-kb  wave,  which  further  indicates 
the  absence  of  rate  effects.  The  resolution  of  this  puzzle  must  await  additional 
measurements  -  preferably  near  the  transition  point. 


SPECIFIC  VOLUME, V - cmVg 

FIG.  5-17  PHASE  DIAGRAM  FOR  IRON  (a  -  hep) 
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(3)  Quartz 


The  two  known  high  density  forms  of  quartz  are  coesite  and  stishovite. 

The  former  is  produced  at  relatively  low  pressures  but  is  a  slow  transition 
and  is  believed  to  be  produced  to  only  a  limited  extent  by  shock  compression 

(28) .  Stishovite  is  believed  to  be  produced  in  shock  (28,29)  under  the  condi¬ 
tions  shown  in  Table  II;  this  is  believed  not  to  be  an  equilibrium  point  (30). 
Nonetheless  the  Wackerle  compression  data  have  been  used  in  Eq.  (5-54)  to 
determine  (dP/dT)cajc  •  The  result,  shown  in  the  last  row  and  column  of 
Table  5 -II  is  about  twelve  times  greater  than  the  static  value  quoted  by  Ahrens 

(29) .  This  can  be  regarded  as  confirmation  of  the  nonequilibrium  character  of 
the  transition.  The  work  by  Ahrens  and  Gregson  suggests  that  shock  recovery 
experiments  can  be  used  to  study  the  kinetics  of  the  transition,  but  this  re¬ 
mains  to  be  done. 

E.  DISCUSSION 

The  ordering  of  slopes  in  the  P-V  plane  for  transitions  of  Types  1  and  2 
appears  to  be  a  useful  minor  tool  in  the  evaluation  of  shock  data.  Transitions 
of  Type  3,  though  curious,  are  of  no  importance.  Of  greater  importance  is 
the  relation  between  adiabatic  slope  in  the  mixed  phase  region  and  dP/dT  ; 
this  allows  cross -comparisons  to  be  made  between  different  kinds  of  measure¬ 
ments.  The  discrepancies  found  between  direct  and  inferred  values  of  dP/dT 
are  so  large  that  needs  for  better  thermodynamic  data  at  high  pressures  and 
for  more  critical  study  of  shock  data  and  their  interpretation  are  evidently 
pressing;  this  clearly  points  the  way  for  further  productive  research. 

The  discrepancies  in  the  iron  data  are  particularly  disturbing.  The 
magnitude  of  the  transition  pressure  agrees  reasonably  well  with  the  static 
value.  The  static  transition  is  reported  to  be  sluggish  (31);  but,  on  the  basis 
of  experience  with  stress  relaxation  effects  (32)  one  would  predict  that  in¬ 
complete  transition  in  shock  would  lead  to  a  first  wave  amplitude  above  the 
static  transition  pressure  and  to  a  decay  of  the  first  wave  amplitude  with  time. 
We  might  assume  that  the  shock  data  are  equilibrium  values  and  that  the 
specific  heat  and  thermal  expansion  coefficients  are  in  error  by  a  large 
amount.  Still  another  possibility  is  that  the  Hugoniot  points  are  at  equilibrium 
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and  beyond  the  mixed  phase  region.  The  discontinuity  at  the  second  phase 
boundary  would  then  lead  to  differences  in  the  observed  direction  between 
^RHl  and  •  Finally  it  is  not  inconceivable  that  the  phase  line  re¬ 

ported  by  Johnson  et  al.  is  seriously  in  error.  It  was  obtained  by  a  clever  but 
unusual  technique  which  is  not  thoroughly  understoc  1;  verification  by  more 
conventional  measurements  would  not  be  out  of  order.  Minshall  made  a 
limited  effort  to  determine  dP/dT  with  pin  techniques,  but  not  all  of  his  data 
have  been  made  available  and  his  conclusions  and  their  certainty  are  not 
known  to  the  authors.  Static  measurements  would  be  important  if  the  pressure 
of  transition  can  be  accurately  determined. 
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6-  FLOW  CALCULATIONS 
J.  O.  Erkman  and  G.  R.  Fowles 


A.  INTRODUCTION 

The  work  on  this  phase  of  the  program  is  intended  to  serve  as  a  guide  for 
assessing  the  sensitivity  of  predictions  of  shock  propagation  from  underground 
explosions  to  uncertainties  in  the  constitutive  relation  (equation  of  state)  of  the 
medium.  In  a  previous  report  on  this  program*  comparisons  were  shown  of  the 
pressure  and  particle  velocity  profiles  predicted  assuming  various  models  for 
the  equation  of  state.  For  these  calculations  spherical  symmetry  was  assumed 
and  the  energy  source  was  taken  to  be  an  adiabatically  expanding  gas  of  pre¬ 
scribed  mass  and  energy.  Variations  of  the  zero-degree  isotherm  (computed 
for  quartz  and  its  high-pressure  polymorph,  stishovite),  GrUneisen's  ratio  and 
the  initial  porosity  were  examined.  Perhaps  surprisingly,  the  differences  in 
the  shock  profiles  showed  no  drastic  effects  of  these  fairly  pronounced  varia¬ 
tions  in  the  equation  of  state.  The  largest  differences  observed  in  the  peak 
pressure  at  a  given  radius,  due  to  combinations  of  parameter  variations, 
amounted  to  factors  of  2  to  3.  Evidently,  some  variations  in  the  equation  of 
state  tend  to  be  self-compensating  in  their  effects  on  shock  propagation.  Tlius, 
in  the  case  of  variable  porosity  the  energy  dissipation  is  known  to  be  highly  de¬ 
pendent  on  porosity,  for  smaller  values  of  porosity,  and  one  might  expect  faster 
decay  of  peak  pressure  in  more  porous  material.  However,  the  high  energy 
loss  per  unit  mass  in  distended  material  compared  with  that  in  initially  compact 
material  tends  to  be  offset  by  the  correspondingly  smaller  total  mass  between 
the  source  and  a  given  radius.  Hence,  the  pressure  is  comparable.  A  similar 
argument  applies  to  variations  in  the  zero-degree  isotherm  in  that  higher  initial 
densities  are  (usually)  associated  with  lower  compressibilities;  again  greater 
total  mass  is  associated  with  less  energy  dissipation  per  unit  mass.  To  these 
effects  must  be  added  the  other  obvious  normalizing  effect,  namely  geometrical 
divergence. 
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Some  of  the  results  obtained  from  the  use  of  the  above  models  are  sum¬ 
marized  in  Figs.  6-10  to  6-13. 

In  view  of  the  previous  results,  an  extensive  and  detailed  study  of  the  pa¬ 
rameters  mentioned  above  was  not  pursued  during  the  past  year.  Instead  effort 
was  concentrated  on  an  area  that  has  so  far  received  little  attention  and  which 
could  cause  peculiarities  in  experimentally  observed  pulse  shapes  and  which 
might  significantly  influence  shock  decay.  This  problem  was  suggested  by  the 
experimental  equation  of  state  measurements  and  is  the  problem  of  the  effects 
of  a  phase  transition  on  shock  propagation.  See  also  Sec.  5.  Polymorphic 
phase  transitions  are  common  for  rocks  and  minerals  because  they  are  gener¬ 
ally  composed  of  open  silicate  structures  that  are  unstable  at  high  presures. 
Moreover,  some  at  least  are  known  to  be  irreversible.  The  high  pressure 
polymorphs  of  quartz ,  coesite  and  stishovite  have  both  been  recovered  in  the 

vicinity  of  underground  nuclear  explosions  and  near  the  Barringer  meteor 
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crater.  Where  a  phase  change  is  irreversible,  large  energy  dissipation  occurs 
that  could  cause  appreciable  attenuation  of  the  shock  pulse.  Perhaps  more 
importantly,  phase  changes  can  cause  the  shock  to  propagate  as  two  distinct 
fronts  of  different  velocity.  The  peak  pressure  can  thereby  be  considerably 
delayed  with  respect  to  the  first  shock  arrival.  Failure  to  recognize  the  exist¬ 
ence  of  such  shock  structure  could  cause  erroneous  interpretation  of  experi¬ 
mental  field  observations ,  such  as  may  be  obtained  with  short  duration  peak 
pressure  gages,  or  simple  time-of-arrival  gages. 

Initial  attempts  to  include  the  phase  transition  in  the  flow  calculations  led 
to  difficulties  in  the  form  of  instabilities.  These  instabilities  are  thought  to 
arise  from  discontinuities  in  slope  of  the  P-V  relation,  although  the  precise 
reasons  for  the  trouble  are  not  understood.  Consequently,  some  effort  was 
devoted  to  generating  functions  that  would  smoothly  represent  both  reversible 
and  irreversible  phase  changes.  The  results  of  flow  calculations  based  on 
these  functions  show  qualitatively  the  kinds  of  pulse  shapes  that  would  be  ex¬ 
pected  in  the  case  of  a  material  that  undergoes  a  phase  transformation.  Even 
with  these  smooth  functions,  however,  oscillations  in  the  flow  parameters  were 
not  entirely  eliminated.  For  the  most  part,  the  phase  change  has  been  con¬ 
sidered  to  be  reversible.  This  case  is  of  interest  because  the  model  implies 
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that  rarefactions  in  the  flow  should  interact  to  form  a  rarefaction  shock.  *  *  * 
Such  a  discontinuity  is  as  difficult  to  treat  with  a  finite  difference  method  as  the 
more  usual  compression  shock.  Hence,  for  a  reversible  phase  change,  it  is 
necessary  to  use  the  artificial  viscosity ,  Q ,  in  the  flow  calculations  for  rare¬ 
factions  as  well  as  for  shocks. 

One  of  the  models  scribed  below  represents  a  reversible  phase  change 
only.  Another  model  w  devised  which  could  be  used  to  represent  both  irre¬ 
versible  and  reversiK  phase  changes.  No  provision  was  made  for  the  phase 
change  of  only  a  portion  of  the  material.  That  is,  the  material  must  be  shocked 
to  a  critical  volume  before  it  is  allowed  to  expand  as  irreversibly.  This  model 
reproduces  the  0°K  isotherm  of  stishovite  at  high  pressure.  Hence,  it  is  rea¬ 
sonable  to  compare  the  results  of  the  flow  calculation  using  this  model  with 
results  reported  previously  for  distended  quartz  and  stishovite. 

B.  EQUATIONS  OF  STATE 

As  noted  in  the  Introduction,  several  relations  have  been  used  for  the  equa¬ 
tions  of  state  in  the  computation  experiments.  In  all  cases,  a  form  of  the  Mie- 
GrUneisen  equation  of  state  has  been  used.  This  equation  relates  pressure  P , 
volume  V,  and  energy  E  for  one  state,  say  on  an  adiabat,  to  a  state  at  the 
same  volume  on  the  0°K  isotherm,  so  that 

P  -  Pk  =  T(E-Ek)/V  (6-1) 

where  subscript  k  refers  to  the  0°K  isotherm  and  T  is  the  Griineisen  ratio. 

In  the  following,  an  equation  of  state  consists  of  a  definition  of  a  0*K  Isotherm, 
which  is  specified  by  a  relation  between  Pk  and  V .  For  porous  materials  the 
initial  state  is  not  a  true  thermodynamic  state  of  the  solid.  For  these  cases, 

therefore,  Pk  is  required  to  be  zero  for  volumes  greater  than  the  initial 

volume  V  .  The  value  of  V  for  quartz  and  stishovite  are  their  respective 
o  o 

crystal  specific  volumes,  0.37  and  0.23  cc/g.  For  dry  playa  (as  reconsti¬ 
tuted  for  the  Hugoniot  experiments  the  value  of  VQ  is  0.513  cc/g. 
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Because  of  the  phase  change  in  quartz,  theoretical  isotherms  were  derived 
for  both  quartz  and  the  high  pressure  material,  stishovite,  see  Fig.  6-1.  The 
phase  change  is  represented  by  a  transition  curve  connecting  the  quartz  to  the 
stishovite  curve  at  a  pressure  of  about  0.2  Mbar.  During  the  first  year's 
work,  *  a  straight  line  was  used  for  the  transition  curve.  This,  and  perhaps 
other  causes,  led  to  instabilities  in  the  computed  results.  This  trouble  made  it 
desirable  to  use  as  simple  a  relation  as  possible  for  the  isotherm  so  that  the 
cause  of  the  instability  could  be  investigated. 

1.  Single  Function  Representation  of  a  Phase  Change 

The  function 

Pk  =  0.2  -  609.  4(V  -  0. 27  )3  +  4114(V  -  0.27)4  (6-2) 


FIG.  6-1  COLD  ISOTHERM  FOR  QUARTZ  AND  STISHOVITE 


is  used  to  relate  the  pressure,  Pk  ,  along  the  0°K  isotherm  to  the  volume  V . 
Units  of  Pk  are  megabars  and  of  volume  are  cc/g.  The  interesting  feature  of 
this  function  is  shown  in  Fig.  6-2.  The  function  is  flat  for  V  =  0.27  cc/g , 
where  the  pressure  is  0.2  Mbar.  The  energy  for  the  0°K  isotherm  is  obtained 
by  evaluating  the  expression 


Ek  ■  •/pkdf  • 

Sound  speed  is  evaluated  by  use  of  the  expression 


(6-3) 


The  initial  volume,  VQ  ,  is  0. 371  cc/g  for  this  model.  Because  the  same  rela¬ 
tion  is  used  for  compressions  as  for  rarefactions,  this  model  represents  a 
reversible  phase  change. 

2.  Multifunction  Representation  of  the  Phase  Change 

A  better  representation  of  the  theoretical  curves  for  quartz  and  stishovite 
can  be  obtained  if  three  functions  are  used.  This  also  permits  the  flow  calcu¬ 
lations  to  proceed  as  if  the  phase  change  were  either  reversible  or  irreversible. 
For  the  interval  0.376  <  V  i  M  where  M  =  0.25,  the  relation 

Pk  =;  0.2  +  (V-0.25)3  [-977.1  +  (V- M)  (12397. 0  -  43129. 0(V - M)| ]  (6-5) 

is  used.  Equation  (6-5)  represents  the  quartz  curve  fairly  well  up  to  about  0. 06 
Mbar,  see  Fig.  6-3.  Above  0.06  Mbar,  the  effect  of  the  phase  change  causes 
Eq.  (6-5)  to  give  results  above  the  quartz  curve.  At  V  *  0,26,  the  value  of 
Pk  is  0. 2  as  desired. 

For  greater  pressures  the  cold  isotherm  is  given  by 
Pk  =  0.2  +  (V-M)3  1-3483.0  +  (V-M)  [-21681.0  -  168022. 0(V-M)|]  (6-6) 
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FIG.  6-2  SINGLE  FUNCTION  REPRESENTATION  OF  PHASE  CHANGE, 
ZERO-DEGREE  ISOTHERM 


SURE  -  Mbar 


FIG.  6*3  MULTIFUNCTION  REPRESENTATION  OF  PHASE  CHANGE, 
ZERO-DEGREE  ISOTHERM 
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This  function  joins  the  stishovite  curve  rather  smoothly  at  V  =  0. 193  ,  see 
Fig.  6-3.  Finally,  for  V  <  0.193  ,  the  stishovite  curve  is  represented  by 


pk  - 


-8.913  +  17  1 374.  23  +7  (5058.0  +  63476  7  I  j 


(6-7) 


where 

n  =  V  -  0.2304. 

Equation  (6-7)  gives  the  same  results  for  pressures  up  to  10  Mbar  as  does  the 
more  complicated  expression  used  previously  for  stishovite. 

The  energy  for  the  cold  isotherm  is  evaluated  by  the  use  of  Eq.  (6-3) ,  using 
the  proper  constants  of  integrations  where  functions  are  joined.  If  the  phase 
change  is  assumed  to  be  reversible,  the  path  ABCD  in  Fig.  6-4  is  used  for  the 
relief  process  as  well  as  for  the  compression  process.  If  the  phase  change  is 
considered  to  be  irreversible,  the  compression  path  is  still  ABCD.  For  the 
relief  process,  however,  the  path  DCE  must  be  followed,  so  that  Ek  =  0  at 
V  =  V2 .  Hence  part  of  the  energy  gained  by  the  material  in  compression  is 
lost  in  the  phase  change. 

It  must  be  understood  that  the  above  scheme  is  not  intended  to  describe  the 
behavior  of  any  real  material.  The  scheme  is  intended  to  be  used  in  flow  cal¬ 
culations  in  order  to  see  what  effects  a  possible  phase  change  has  on  the  com¬ 
puted  results.  This  model  makes  no  provision  for  a  change  in  phase  of  a  part 
of  the  material.  In  order  to  expand  as  stishovite,  the  material  must  have  been 
compressed  to  a  volume  of  0. 193  cc/g  or  less,  for  which  the  corresponding 
pressure  is  0,72  Mbar. 

C.  SPHERICAL  PISTON 

The  code  simulates  a  situation  in  which  a  disturbance  in  the  playa  is  driven 

by  an  event  in  a  cavity  1  meter  in  radius.  Thus  the  cavity  has  about  the  same 
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volume  as  the  cavity  for  the  Rainier  event.  The  cavity  is  assumed  to  be  filled 
with  an  ideal  gas  whose  density  is  0. 14  g/cc  and  whose  energy  is  73. 5  megabars 
cc/g  (or  73.5  x  10  ergs/g)  so  that  the  pressure  is  about  6.9  megabars.  The 
adiabatic  exponent  is  assumed  to  be  5/3.  The  advantage  of  using  this  model  is 
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FIG.  6-4  SCHEMATIC  OF  REVERSIBLE  AND  IRREVERSIBLE 
PHASE  CHANGE 


that  the  maximum  pressure  induced  in  a  medium  surrounding  the  cavity  depends 
on  the  P,V,E  relation  of  the  media.  In  some  of  the  earlier  work,  an  arbitrary 
relation  was  assumed  between  the  pressure  in  the  cavity  and  the  time,  so  that 
the  same  maximum  pressure  was  induced  in  a  "soft"  material  as  in  a  "hard" 
material. 

The  assumption  that  the  bomb  fills  the  cavity  uniformly  with  an  ideal  gas 
which  expands  adiabatically  precludes  the  formation  and  interaction  of  shock 
waves  in  the  cavity.  In  the  calculations  the  cavity  is  not  zoned,  and  no  details 
of  the  flow  of  the  gas  in  the  cavity  are  calculated. 
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D.  RESULTS  OF  CALCULATIONS 

Results  of  flow  calculations  will  be  presented  here  for  each  of  the  equations 
of  states  described  in  Section  B.  Results  are  presented  graphically,  largely  by 
showing  the  dependence  of  both  the  pressure  and  the  particle  velocity  on  distance 
at  fixed  times  for  a  particular  equation  of  state.  Time  of  arrival  curves  and 
plots  of  peak  pressure  vs  distance  are  also  used  to  show  some  of  the  results  of 
the  changes  in  the  equations  of  state. 

1 .  Results  Using  Single  Function  Representation  of  a  Phase  Change 

Equation  (6-2)  is  combined  with  the  Mie-Griineisen  equation  of  state,  Eq. 
(6-1),  the  value  of  T  being  2/3.  Calculations  were  performed  for  initial  vol¬ 
umes  of  0. 371  cc/g  and  for  0. 504  cc/g.  For  these  calculations,  the  bomb  energy 
was  reduced  from  73.  5  to  7. 35  Mbar  cc/g  so  that  the  maximum  pressure  in  the 
cavity  is  0. 7  Mbar. 

Results  obtained  from  the  flow  calculations  for  an  initial  volume  of  0. 371 
cc/g  on  are  presented  in  Figs.  6-5(a)  and  6-5(b).  Figure  6-5(a)  shows  the 
pressure  as  a  function  of  distance  from  the  center  of  the  cavity.  The  interesting 
features  of  the  pressure  profiles  are  that  the  fronts  of  many  of  them  are  no 
steeper  than  portions  of  the  back  sides  of  the  waves.  For  these  calculations, 
the  artificial  viscosity  term  was  employed  in  both  compression  and  in  rarefac¬ 
tion  waves.  This  is  consistent  with  the  original  formulation  of  the  artificial 
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viscosity  method  (except  that  a  small  linear  term  has  been  introduced). 

Further  complications  result  when  the  material  is  distended,  i.e. ,  the 
quartz  is  assumed  to  be  filled  with  small  holes.  It  is  also  assumed  that  when 
the  temperature  is  0°K,  the  materia)  crushes  under  little  or  no  pressure. 

Figures  6-6(a)  and  6-6(b)  represent  the  results  of  computations  for  which  the 
original  volume  is  0. 50  cc/g  or  about  1. 36  times  the  volume  used  in  obtaining 
the  results  discussed  immediately  above.  The  results  are  very  different  in  that 
the  pressure  drop  in  the  initial  rarefaction  is  very  steep.  When  the  time  is 
200  gsec,  the  rarefaction  has  almost  overtaken  the  shock  front.  When  the  time 
is  250  gsec,  the  entire  wave  has  been  drastically  attenuated.  Following  this 
initial  attenuation,  the  wave  is  oscillatory  back  of  the  shock  front.  From  the 
calculations  themselves,  there  is  no  way  of  determining  if  the  oscillatory 
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FIG.  6-5(b)  PARTICLE  VELOCITY  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 

Single  function  phase  change,  Vq  -  0.37  cc/g;  \R  5  cm,  Q  operated  in  compression  and  rarefaction. 


FIG.  6-6(a)  PRESSURE  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 

Single  function  phase  change,  VQ  =  0.50  cc/g,  Ar  =  2.5  cm,  Q  operated  in  compression  only. 


0250 


FIG.  6-6(b)  PARTICLE  VELOCITY  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 

Single  function  phase  change,  Vg  =  0.50  cc/g,  Ar  =  2.5  cm,  Q  operated  in  compression  only, 


behavior  is  due  to  wave  interaction,  or  if  it  results  from  deficiencies  in  the 
method  of  obtaining  the  solution.  The  undershoot  on  the  early  profiles  (T  = 
100,  150  and  200  p sec)  is  likewise  not  understood. 

The  results  presented  in  Figs.  6-6(a)  and  6-6(b)  were  obtained  by  using 
Q  in  the  conventional  manner,  i.e. ,  in  the  "cut-off1  form  so  that  Q  is  non¬ 
zero  only  in  compression.  If  Q  is  used  as  originally  proposed  by  von  Neumann 
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and  Richtmyer  the  results  shown  in  Figs.  6 -7 (a)  and  6 -7(b)  are  obtained.  For 
these  results,  there  are  no  oscillations  following  the  initial  large  rarefaction. 
Neither  are  there  undershoots  of  pressure  or  of  particle  velocity  immediately 
following  this  rarefaction.  The  profiles  are  still  oscillatory,  however,  after 
the  pressure  has  been  attenuated  below  100  kbar.  Hence  the  significance  of 
these  oscillations  is  still  undetermined. 

The  use  of  a  linear  artificial  viscosity  should  have  a  significant  damping 
effect  on  oscillations.  The  term 

Q  =  -0.5(u(J+l)  -  u(J)) 

was  used  alone,  and  only  in  compression  in  obtaining  the  results  shown  in  Figs. 
6-8(a)  and  6-8(B).  These  results  have  ail  the  peculiarities  first  noted  in  Figs. 
6-6(a)  and  6-6(b).  That  is,  the  solution  undershoots  following  the  first  large 
rarefaction  and  then  oscillates.  Also,  after  attenuation  to  below  100  kbar  the 
profiles  are  oscillatory. 

From  the  three  sets  of  results  as  shown  in  Figs.  6-5(a)  through  6-8(b),  it 
appears  that  interactions  between  the  form  of  the  equation  of  state  and  the  nqmer 
ical  computation  scheme  may  be  causing  serious  convergence  problems.  The 
oscillation  and  undershooting  mentioned  above  are  also  observed  when  the  bomb 
energy  is  increased  to  73. 5  Mbar  cc/g.  Thus  these  features  of  the  calculated 

results  are  not  dependent  on  the  pressure  in  the  cavity.  If  any  judgment  is  to  be 
made,  the  results  shown  in  Fig.  6-7  are  preferable.  That  is,  permitting  the 
artificial  viscosity  to  operate  throughout  compressions  and  rarefactions  give 
the  preferred  results. 


FIG  6-7(a)  PRESSURE  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 

Single  function  phase  change,  VQ  =  0.50  cc/g,  AR  =  2.5  cm,  Q  operated  in  both  compression  and  rarefaction 


FIG.  6-7(b)  PARTICLE  VELOCITY  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 

Single  function  phase  change,  Vq  =  0.50  cc/g,  AR  =  2.5  cm,  Q  operated  in  both  compression  and  rarefaction. 


FIG.  6-8(b)  PARTICLE  VELOCITY  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 

Single  function  phase  change,  Vn  =  0.50  cc/g,  AR  =  2.5  cm,  Q  =  —0.5  [u(JH)  —  u(J)]  in  compression  only. 


2.  Multifunction  Representation  of  a  Phase  Change 

As  described  in  Section  B,  this  representation  of  the  ,  V  relation  rep¬ 
resents  quartz  at  low  pressure  and  stishovite  at  high  pressure.  In  this  formu¬ 
lation  the  phase  change  can  be  treated  as  if  it  were  completely  reversible ,  or 
stishovite  may  be  allowed  to  expand  entirely  as  stishovite,  but  no  mixed  phases 
are  possible  in  expansion.  The  critical  compression  for  which  alternate  relief 
paths  may  be  chosen  corresponds  to  V  =  0.193  cc/g. 

These  equations  of  state  relations  have  been  used  in  calculations  for  which 
the  original  volume  is  0.  51  cc/g  (that  is,  the  material  is  distended  from  the 
solid  specific  volume  of  0.371  cc/g),  and  for  which  the  cell  size  is  5  cm.  The 
energy  in  the  cavity  is  73. 5  Mbar  cc/g,  giving  a  maximum  pressure  of  about  7 
Mbar  in  the  material  near  the  cavity.  The  artificial  viscosity  term  is  used  in 
both  compression  and  in  rarefaction.  Results  are  shown  in  Figs.  6-9(a)  and 
6-9(b)  for  which  all  elements  of  the  compressed  i.iaterial  expand  reversibly 
with  respect  to  the  cold  isotherm.  This  model  also  shows  the  steep  rarefaction 
(see  profiles  for  200  and  300  jisec)  discussed  above.  The  particle  velocity  pro¬ 
files  are  perhaps  more  interesting  for  this  case  than  are  the  pressure  profiles. 
As  stated  above,  there  is  no  intent  to  imply  that  playa  (or  any  other  material)  is 
described  by  the  models  discussed  here.  However,  the  response  of  some  rea 
material  may  have  features  in  common  with  a  model  such  as  the  one  used  here. 

The  results  obtained  when  stishovite  is  require.!  to  expand  entirely  as 
stishovite  differ  very  little  from  the  results  shown  in  Fig.  6-9.  For  this 
reason,  the  profiles  are  not  given.  The  explanation  of  this  small  difference  is 
that  very  little  of  the  material  is  ever  compressed  to  a  volume  of  0. 193  cc/g 
or  less.  Out  to  about  150  cm  from  the  center  of  the  cavity,  some  of  the  mate¬ 
rial  is  compressed  so  that  the  volume  is  less  than  0. 193  cc/g.  Slightly  further 
out  the  wave  has  attenuated  to  such  an  extent  that  the  volume  is  never  reduced 
to  the  critical  volume.  Thus  a  spherical  shell  50  cm  thick  surrounding  the  100- 
cm  radius  cavity  contains  stishovite.  Forcing  this  small  amount  of  material  to 
expand  as  stishovite,  changes  the  propagation  of  the  wave  very  little.  Using  a 
greater  pressure  in  the  cavity  would  undoubtedly  make  the  results  obtained 
from  the  use  of  the  two  models  differ  more  dramatically  in  range  of  distances 
examined. 


FIG.  6-9(a)  PRESSURE  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 

Multifunction  phase  change  (reversible)  Vq  =  0.51  cc/g,  AR  =  5  cm, 

Q  in  compression  and  rarefaction. 

E.  CONCLUSIONS 

In  this  phase  of  the  program,  including  the  previous  years  effort,  some 
effects  due  to  variations  in  equation  of  state  on  spherical  shock  propagation  have 
been  examined.  For  this  study  a  particular  energy  source  was  chosen  to  give  a 
rough  approximation  to  a  representative  nuclear  explosion.  Pressure  and  par¬ 
ticle  velocity  pulse  shapes  and  the  rates  of  decay  of  their  peak  amplitudes  with 
distance  are  shown  for  different  models  of  the  equation  of  state  based  on  differ¬ 
ent  assumptions  for  the  zero-degree  isotherm,  Griineisen's  ratio,  the  degree  of 
porosity,  and  including  a  reversible  phase  change.  These  variations  represent 
reasonable  bounds  to  the  uncertainties  in  knowledge  of  the  equation  of  state  of 
playa.  Because  of  the  particular  source  function  chosen,  no  claim  of  generality 


135 


FIG.  6-9(b)  PARTICLE  VELOCITY  vs.  DISTANCE  FROM  CENTER  OF  CAVITY 
Multifunction  phase  change  (reversible)  Vq  =  0.51  cc/g,  AR  =  5  cm, 

Q  in  compression  and  rarefaction. 


of  the  results  can  be  made.  Nevertheless,  the  results  should  be  useful  in  indi¬ 
cating  the  sensitivity  of  shock  propagation  to  uncertainties  in  the  equation  of 
state. 

The  effects  of  changes  in  chemical  composition  within  the  limits  of  observed 
variations  in  the  playa  were  not  studied  explicitly  by  shock  calculations.  Exam¬ 
ination  of  the  effects  on  the  equation  of  state  showed  that  the  variation  (theoreti¬ 
cally)  is  small  and  lies  within  the  range  of  variation  already  examined  by  means 
of  other  parameters,  such  as  Griineisen’r  ratio.  It  is  concluded,  therefore, 
that  variations  in  chemical  composition  are  of  lesser  importance  than  the  uncer¬ 
tainties  to  be  expected  in  other  parameters. 
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Unfortunately,  experimental  data  on  the  effects  of  moisture  on  the  equation 
of  state  (particularly  release  adiabats)  did  not  become  available  in  time  to  in¬ 
clude  such  cases  in  the  shock  calculations.  Moreover,  a  strictly  theoretical  ap¬ 
proach  did  not  appear  likely  to  be  very  meaningful.  Consequently,  the  effects  of 
moisture  on  shock  propagation  remains  an  unstudied  problem. 

The  observed  differences  in  the  arrival  time  of  the  peak  pressure  are  shown 
in  Figs.  6-10  and  6-11.  The  differences  in  peak  pressure  as  function  of  the 
radial  distance  amount  to  a  factor  of  three  or  iess,  and  are  shown  in  Figs.  6-12 
and  3-13.  In  general  the  peak  pressures  assuming  a  cold  isotherm  correspond¬ 
ing  to  stishovite  lie  below  those  assuming  a  quartz  isotherm.  The  direction  of 
this  result  is  clearly  proper  since  there  is  greater  energy  dissipation  (waste 
heat)  where  the  release  adiabats  (and  cold  isotherm)  are  steeper. 

The  effects  of  varying  Griineisen's  ratio  are  also  in  the  expected  direction. 
Higher  values  of  Griineisen’s  ratio  correspond  to  shallower  adiabats  and,  hence, 
less  energy  dissipation. 

The  effects  of  varying  porosity  are  surprisingly  small  because  porosity 
drastically  effects  energy  dissipation.  The  only  plausible  explanation  is  that 
there  are  two  effects  that  tend  to  be  self-compensating.  Higher  porosity,  and 
higher  dissipation  arc  associated  with  less  mass  between  the  source  and  a  given 
radius.  These  two  effects  cancel  to  such  a  degree  that  the  peak  pressure  at  a 
given  radius  is  not  seriously  affected.  This  result  is  perhaps  the  most  signifi¬ 
cant  one  of  this  phase  of  the  program. 

The  effect  of  a  phase  change  is  mostly  to  alter  the  pulse  shape  and  to  cause 
abrupt  changes  in  the  rate  of  decay  of  peak  pressure.  Recognition  of  such  pos¬ 
sibilities  may  be  important  in  understanding  field  results  obtained  by  short 
duration  pressure  gages,  or  simple  time-of-arrival  gages. 

It  was  also  found  that  where  reversible  phase  changes  occur,  the  artificial 
viscosity  term  in  the  numerical  differencing  scheme  should  be  allowed  to  oper¬ 
ate  in  rarefaction  as  well  as  compression  in  order  to  accommodate  rarefaction 
shocks. 
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FIG.  6-12  PEAK  PRESSURE  AS  FUNCTION  OF  RADIUS  BASED  ON  QUARTZ  EQUATION  OF  STATE 
(See  Fig.  6~10  for  definition  of  symbols) 
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HIGH  DENSITY  MODELS  FOR 

SEMIEMPIRIC AL  COLD  EQUATIONS  OF  STATE  FOR  SOLIDS 


Christian  Peltzer 


I.  INTRODUCTION 

Various  high  density  models  for  solids  have  been  introduced  to 
calculate  equations  of  state  of  materials  valid  in  the  pressure  range 
attainable  through  explosions,  high  velocity  impacts,  etc.  It  is  impor¬ 
tant  to  realize  that  the  resulting  equations  of  state  are  "asymptotic 
model  equations  of  state,  "  not  proper  mathematical  asymptotic  forms 
(in  the  limit  of  very  high  densities)  of  the  exact  equation  of  state  of  the 
system.  Because  of  their  ad  hoc  nature  it  is  difficult,  if  at  all  possible, 
to  establish  their  range  of  validity  and  to  give  estimates  of  the  error 
involved  in  their  use  in  various  density  ranges.  An  answer  to  such 
questions  can  only  be  obtained  within  a  more  general  framework  that 
would  allow  a  precise  mathematical  formulation  of  the  physical 
assumptions  made  and  a  bona  fide  estimate  of  the  resulting  error. 

The  basic  physical  assumptions  underlying  all  high  density  model 
are  that,  in  the  limit  of  very  high  densities,  all  shell  structure  of  in¬ 
dividual  atoms  disappears  and  that  the  energy  of  the  electron  system 
tends  towards  that  of  a  uniform  density  system. 

Here  we  discuss  briefly  various  high  density  models*  and  the 
corresponding  cold  equation  of  state. 


*  Only  nonrelativistic  models  will  be  considered. 


A  general  and  powerful  method  for  handling  the  basic,  many -body 
problem  of  solid  state  physics  has  been  developed  by  several  authors. 

An  outline  of  this  method,  the  generalized  density  operator  theory,  will 
be  given  elsewhere;  there,  we  will  try  to  show  how  it  allows  a  unification 
of  various  cohesive  energy  theories  and  a  better  understanding  of  the 
approximations  involved,  and  why  it  could  lead  to  new  and  better  compu¬ 
tational  techniques  for  obtaining  equations  of  state.  Any  serious  discussion 
of  the  validity  of  the  various  models  used  below  is  best  deferred  until  this 
formalism  has  been  introduced. 

II.  SIMPLEST  HIGH  DENSITY  MODELS 
A  .  Perfect  Free -Electron  Gas 

The  simplest  high  density  model  (Model  I)  is  that  of  a 
perfect,  free  (i.e.  ,  noninteracting),  electron  gas  of  uniform  density  in 
a  uniform  positive  background  (necessary  for  charge  neutrality).  Such 
a  system  possesses  only  kinetic  energy,  the  constant  potential  being 
taken  as  zero;  this  energy  is  a  function  of  the  electron  density  and  so, 
for  a  fixed  number  of  electrons,  of  the  volume.  (The  electrons  being 
fermions  cannot  all  lie  in  the  same  lowest  one-electron  energy  level 
even  at  T  =  0 . ) 


*  Forthcoming  ISR  report  ( 18453 1  -  1 10) . 
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Let  V  be  the  specific  volume  of  the  system,  M  the  molecular 

m 

0 

weight,  Z  =  Sn.Z.  the  total  number  of  electrons  per  molecule,  p  the 
electron  density.  Then 


N  Z 

e  a 


MV  MV 

m  m  oo 


N  z  7  .1 

£  -  — —  £  =  const  (1) 


where  N  is  Avogadro's  number,  £  =  -  the  relative  specific  volume 

a  V 

oo 

(V  =  V  at  the  reference  state  P  =  0,  T  =  0),  and  v  the  effective 
oo  a 

atomic  volume  at  the  reference  state. 


The  electron  kinetic  energy  density  £  of  a  perfect  free  electron 

K 


gas  is  well  known  and  is 


5/3 


fk =  \p 


(2) 


with 


Oi  = 
k 


3  Ww2 

IT  |3,r 


2/3  ?  -27  2 

—  =  2.87le  a  =  3.  505  x  10  erg  cm  (3) 
m  o 


n  -8 

where  a  =  - —  =  Bohr  radius  =  0.  5292  x  10  cm. 

o  2 

me 


The  internal  energy  per  gram  E  =V£  is  given  by 

C  K 


with 


<■  v‘n-v“’ 


A.  =  a. 
k  k 


/  N  Z  \  5/3 
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=  a. 
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5/3 
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-2/3 

B,  =  V  A .  =  a,  V 
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(5) 


The  corresponding  cold  pressure  and  bulk  modulus  are  obtained  from 
the  the  rmodynamic  relations: 
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(6) 


P 
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and  so 


is  obtained. 


dE 


dE 


dV 


V  dt 
oo 


dP  dP 

BC  ■  '  v"dv£  =  '  ?  SC 


p1  =  4  a,v-5/3 

c  3  k 


■  fv-1  b,  r5/3 

3  oo  k 


(7) 


b1  =  4Akv5/3  =  4v‘  B.rs/3 

c  9  k  9  oo  k 


The  following  relations  should  be  noted 

PIv  =  i.  E1 

c  3  c 


(Virial  theorem) 


B1  ■  *  P1  . 

c  3  c 

This  model  obviously  does  not  exhibit  any  binding  as  E*  is  a  monotonic 
increasing  function  of  p  =  1/V  in  the  whole  domain  (0,°°).  The  energy 
zero  chosen  here  corresponds  to  an  infinitely  dilute  system,  i.  e.  ,  to 
the  limit  p  — ►  0  or  V  — . 


B.  Uniform  Electron  Gas  with  Exchange 

Apparently,  a  better  model  could  be  that  of  a  uniform 
density,  imperfect,  free  electron  gas  in  a  uniform  positive  background. 
By  imperfect  we  mean  that  the  electrons  are  allowed  to  interact  among 
themselves;  the  system  will  then  have,  besides  its  kinetic  energy,  a 
volume  dependent  potential  energy. 


The  total  energy  density  can  be  written  as  a  sum  of  three  terms 


C  =  C,  +  €  +  € 

k  ex  corr 


(8) 


where  the  first  term  is  the  kinetic  energy  of  the  corresponding  free 

electron  gas  (also  often  written  C  and  called  the  mean  Fermi 

r 
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energy*),  the  second  term  is  the  exchange  energy, and  the  third  one  the 
correlation  energy. 

The  exchange  effect  is  of  a  repulsive  nature;  its  origin  lies  in  the 
Pauli  exclusion  principle  that  prevents  two  electrons  with  parallel 
spins  from  occupying  the  same  cell  in  phase  space.  (The  Fermi  statis¬ 
tics  only  limit  the  number  of  electrons  per  cell  to  two  regardless  of 
their  spins.  )  It  follows  that  an  electron  effectively  will  strongly  repel 
all  other  electrons  having  parallel  spins,  and  one  says  that  each  electron 
is  accompanied  by  a  "Fermi  hole.  " 

The  exchange  energy  of  a  free  electron  gas  has  often  been  calcu¬ 
lated  and  is  given  by: 


€  =  -  a  p 

ex  x 


with 


O  3  \ 1 /3  2  2  -19 

0!  =  —  —  I  e  =  0.7386  e  =  1.704  x10  erg  cm  . 

Neglecting  the  Coulomb  correlation  effects,  the  equation  of  state  of 
Model  II  is  given  by 

■  v'-2'3.v-1'3 .  Bkr2/3-Bxr1/3  no, 


with 


N  Z 

A  =  a  r~ 
x  x  M 
i  m 


x  l  v 


B  =  v  -1/3a  =  a  v  — 

X  oo  X  X  OO  V 


*  One  must  not  confuse  the  energy  at  the  Fermi  level  €  and  the 

3  F 

mean  Fermi  energy  c  ,  here  equal  to  —  (  The  same  remark 

F  j  F 

holds  for  the  correlation  energies. 
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From  Eqs.  (10)  and  (6)  we  obtain 


and 
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It  is  easily  verified  that  the  Virial  theorem  also  holds  for  this  model. 


p^V  =  -L  E  +  _L  E 

c  3  kin  3  pot 


where  E.  .  -  E^  and  E  =  -E  =  exchange  energv- 

kin  c  pot  x  6  67 


Model  II  corresponds  to  the  Hartree-Fock  approximation  for  a 
uniform  density  electron  gas;  the  kinetic  energy  of  the  system  is  the 
same  as  that  of  a  perfect  electron  gas,  but  this  model  exhibits  some 
binding  because  of  the  repulsive  nature  of  the  exchange  interaction 
which  is  included  now.  Its  equilibrium  volume  ,  i.e.  ,  the  volume 


minimizing  E^,  is  obtained  as  the  solution  of 

P11  =  0  , 

C 


equ 


and  one  has 
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.  95 
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24 


or 


(13) 


=  41. 95 - — 

’equ  V  M 


=  69.622 


oo  m 


v  x  10 
a 


24 


II 


One  sees  that  £  “  is,  in  general,  much  larger  than  1;  the  binding  is 
equ 

too  weak.  At  the  reference  state,  Model  II  gives  a  pressure  of  the 
order  of  the  megabar  instead  of  zero. 
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Using  Eqs  (10)  and  (12)  as  semiempiric  al  equations  of  state, 
treating  A as  an  adjustable  parameter  to  be  determined  from  the 


relation 


II 

equ 


1 


2A, 


V  i  A 

OO  \  X 


=  1, 


is  essentially  McMillan's  empirical  procedure  applied  to  the  uniform 
electron  gas  model  instead  of  the  T-F  model. 

C .  Uniform  Electron  Gas  with  Exchange  and  Correlation 

The  origin  of  the  correlation  energy  lies  in  the  Coulomb 

interaction  of  the  electrons;  this  is  a  classical  two-body  interaction 

2 

given  by  terms  of  the  form  e  /r..  where  r..  is  the  distance  between 


,th 


th 


U 


U 


the  i  and  j  electrons.  It  depends  on  the  simultaneous  positions  of 
the  two  electrons,  hence  the  name  "correlation  effects."  The  Coulomb 
interaction  between  like  sign  charges  being  of  a  repulsive  nature  and 
this  repulsive  energy  becoming  infinite  for  r  =  0,  one  can  say  that 
each  electron  is  surrounded  by  a  "Coulomb  hole"  with  respect  to  a*l 
other  electrons.  As  the  exchange  effect  gives  rise  to  a  similar  effect 
for  each  electron  with  respect  to  all  other  electrons  with  parallel  spins, 
the  main  part  of  the  Coulomb  correlation  energy  of  an  electron  gas  will 
arise  from  the  interaction  between  electrons  with  opposite  spins  if  the 
exchange  effect  has  already  been  taken  into  account.  The  Coulomb 
correlations  give  rise  to  a  change  in  both  the  kinetic  energy  of  the 
electrons  (with  respect  to  that  of  a  perfect  free  electron  gas)  and  their 
potential  energy,  the  net  effect  being  a  lowering  of  the  total  energy  of 
the  system,  i.  e.  ,  an  increase  in  the  binding. 

A  clear  distinction  must  be  kept  between  a  uniform  density  elec¬ 
tron  gas  and  a  uniform  continuous  charge  distribution.  If  we  replace 
the  electron  gas  by  such  an  equivalent  (charge  wise)  continuous  distri¬ 
bution,  the  electrostatic  potential  energy  of  the  total  system  (electron 
and  nuclei)would  be  zero  here,  each  volume  element  being  neutral. 
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The  electrostatic  potential  energy  of  the  continuous  charge  distribution 
equivalent  to  the  electrons  can  be  called  the  "classical  electrostatic 
potential  energy"  of  the  electrons,  But  in  such  a  replacement  the 
localized  nature  of  the  electrons  is  neglected,  and  the  direct  electron- 
electron  interaction  is  only  partly  taken  into  account,  appearing  as  a 
classical  screening  effect.  The  (Coulomb)  correlation  energy  will  be 
defined  here*  as  the  difference  between  the  energy  calculated  on  the 
basis  of  the  Hartree-Fock  approximation  and  the  exact  one  for  a  given 
Hamiltonian. 


The  evaluation  of  the  average  Coulomb  correlation  energy  density 

T  is  a  considerably  more  difficult  problem  than  that  of  the  exchange 

corr 

energy,  and  explicit  expressions  valid  for  all  densities  have  not  yet 

been  obtained.  Wigner  derived  expressions  for  T  in  the  low  den- 

cor  r 

sity  limit,  and  more  recently  Macke,  Bohm  and  Pines,  and  others 
have  given  high  density  limit  expressions.  Various  interpolation  for¬ 
mulas  have  also  been  used  by  several  authors.  Here  we  shall  only 
consider  the  high  density  limit  expression  in  view  of  the  asymptotic 
nature  of  the  model  and  its  general  shortcomings. 


One  has  for  c  : 

c 


_  e .  .  e, 1 /3  e 

c  =  -  a  p  in  a(p  )  -  a  ,p 

c  c  d 
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a  =  0.853x  10  cm 
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a,  =  0.048  — —  =  2.092x  10  erg 
d  a 
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*  Several  other  definitions  of  the  correlation  energy  are  used  in  the 
lite  rature . 
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Model  III  is  that  of  a  uniform  density,  imperfect  electron  gas  in 
a  uniform  positive  background,  exchange  and  correlation  effects  being 
taken  into  account  in  an  approximate  way  (essentially  through  a  high 
density  perturbation  expansion). 

One  obtains  for  the  cold  internal  energy  per  gram  : 

E111  =  A  V'2/3  -  A  V"1/3  +  A  lnV  -  A 
c  k  x  c  d 


or 


(15) 
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c  k  x  c  d 
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The  cold  pressure  and  cold  bulk  modulus  are  given  here  by 


P111  =  -5/3  1  ^  y-4/3 
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The  inclusion  of  the  correlation  increases  the  binding,  but  one  can 
verify  that  the  pressure  at  the  reference  state  is  still  of  the  order  of  the 
megabar  in  general.  The  equilibrium  volume  is  given  by  the 


positive  root  of  the  equation  P 


Ml 


III 


equ 


equ 


=  0,  i.  e.  , 
4A. 


1  3 


L  A  + 
x 


V 


A  +  8A.  A  J 
x  k  c 


(19) 


and,  in  general,  the  corresponding  is  >>  1  . 


Because  of  the  (Xn  V)-term,  Eq  (17)  cannot  be  used  over  the  whole 
range  (O,00)  of  V.  To  obtain  an  expression  valid  in  this  whole  range  one 
could  simply  replace  the  €  ^  of  Eq  (14)  by  an  interpolation  formula 
having  (14)  as  high  density  asymptotic  limit  and  going  over  into  Wigner's 
expression,  e.g.,  in  the  limit  of  low  densities.  The  validity  of  such  an 
interpolation  formula  is,  of  course,  always  open  to  question  at  inter¬ 
mediate  densities.  (See  II). 

13,14 

D.  Uniform  Electron  Gas  with  a  Positive  Lattice  Background 


In  the  previous  models,  charge  neutrality  of  the  system  was 
obtained  by  assuming  the  existence  of  a  uniform,  continuous,  positive 
charge  distribution;  such  a  background  differs  markedly  from  that  of 
real  matter,  and  it  seems  more  natural  to  represent  the  nuclei  system 
by  a  positive  point  charge  distribution.  The  model  to  be  considered  now 
is  that  of  a  uniform  electron  gas  with  such  a  positive  lattice  background. 

The  essential  assumptions  made  are  that  the  kinetic  energy  of  the 
nuclei  is  negligible,  that  the  static  approximation  is  valid  (allowing  the 
separation  of  the  electrons  and  nuclei  motions),  and  that  the  nuclei 
form  a  regular  lattice.  For  simplicity  we  will  consider  here  only  one 
kind  of  nuclei  and  a  simple  lattice  (no  basis). 
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The  assumption  of  a  regular  lattice  is  made  plausible  (for  all 

elements  except  He  and  possibly  Li)  by  the  fact  that  the  energy  associated 

14 

with  the  formation  of  the  lattice  is  negative  and  by  Abrikosov's  estimate 
of  the  zero  point  nuclear  vibration  amplitudes  for  this  model.  As  Z, 
the  atomic  number,  factors  out  in  the  expression  for  the  lattice  binding 
energy,  the  geometry  of  the  energetically  most  favorable  lattice  is  in¬ 
dependent  of  Z  (universal  for  all  elements).  In  principle  it  should  be 
obtained  by  strict  maximization  of  this  lattice  energy,  but  the  energy 
differences  for  different  lattices  are  small  and  within  the  limits  of 
accuracy  of  the  model  one  can  limit  oneself  to  a  cubic  one.  On  the 
other  hand,  for  polyatomic  materials  energy  differences  between  vari¬ 
ous  nuclei  geometries  could  be  greater,  and  further  study  of  this  point 
seems  worthwhile. 

The  total  energy  of  the  system  can  be  written  as  a  sum  of  three 
terms,  and  in  the  high  density  limit  the  leading  one  is  the  kinetic  energy 
of  the  corresponding  perfect  free  electron  gas  €  .  In  the  approximation 

K 

considered  the  exchange  and  correlation  energies  will  be  the  same  as 
before,  up  to  an  additional  electron-nuclei  correlation  energy.  But  the 
"classical"  electrostatic  energy  of  the  electrons,  € ^  ,  will  no  longer 
be  zero  (as  with  a  uniform  positive  background).  There  are  three  con¬ 
tributions  to  €  the  nuclei-nuclei,  the  electron-nuclei,  and  the 
cl 

electron-electron  Coulomb  interaction  energies;  but  the  three  must  be 
handled  together  as  each  separately  would  diverge  for  an  infinite  solid. 

We  emphasize  again  that  in  calculating  €  ^  the  electrons  are  replaced 
by  continuous  charge  distribution  of  density  (-ep  )  and  that  the  electron- 
electron  interaction  taken  into  account  is  only  the  classical  screening 
effect  of  such  a  charge  distribution,  i.e.,  onlv  a  partial,  average, 
electron-electron  Coulomb  interaction. 
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One  obtains*  for  €  ,  : 

cl 


>4/3 


fcl  =  -  “clP 


with 


a  ^  1.44e2ZZ/3.  (20) 


This  energy  results  from  the  formation  of  the  lattice  and  is  sometimes 
called  the  "lattice  binding  energy.  " 


Neglecting  correlation  terms  we  obtain  the  Model  IV  equation  of 


state 


or 


EIV  s  A  v'2/3-(A  +  A  ,)V1/3 
c  k  x  cl 


EIV  =  B.r2/3  -  (B  +B  l)C-1/3 


(21) 


with  A,  ,  A  ,  B,  ,  B  given  as  before,  and 
k  x  k  x 


A  ,  -  1.44eV/3 

cl 


and 


B  =  A  V  "1/3  =  1.44e2V  Z2/3 

cl  cl  oo  oo 


N  z  i 

4/3  ?  , yi 

1  ZV 

a 

=  1.44e2Z2/3 

oo 

M  j 

1  v 

m  1 

i  a 

4/3 


(22) 


4/3 


Comparing  (21)  and  (10)  (uniform  electron  gas  with  exchange  in  uniform 
positive  background),  one  sees  that; 

1) 


2) 


The  lattice  binding  energy  E  is  of  the  same  order 

in  V  as  the  exchange  energy  (in  the  high  density  limit). 

2 

E  is  proportional  to  Z  while  E  is  proportional 

cl  4/3  x 

to  Z  only;  one  has 


Ex  Y  z-2/3_  0_74  z-2/3^514z-2/3 


c  1 


a 


cl 


1. 44 


(23) 


14  . 


*  Abrikosov  gives  =  1.3;  we  obtained  1.41  for  a  sc  lattice  and 
1.44  for  a  fee  and  a  bcc. 
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3) 


We  can  also  consider  a  Model  V  by  introducing  here  the  correlation 
energy;  these  corrective  terms  are  of  the  same  form,  in  the  high  den¬ 
sity  limit,  as  in  the  case  of  the  uniform  background,  and  their  effect 
is  similar  (see  ID- 


As  we  have  seen,  Models  II  and  III  exhibit  considerably  too  weak 
a  binding,  while  Model  IV  exhibits  too  strong  a  one  at  normal  densities. 
Although  these  two  models  could  be  used  as  bounds  for  the  exact  cold 

isotherm  E  ,  such  a  property  has  in  no  way  been  proved  here.  The 

c 

rigorous  derivation  of  sharp  enough  upper  and  lower  bounds  for  E^,  P^, 

B  over  an  extended  range  of  densities  appears  to  be  an  almost  impos- 
c 

sible  problem  in  the  frame  of  the  standard  approaches  to  the  study  of 
the  ground  state  of  most  quantum  mechanical  systems;  but  we  believe 
that  significant  progress  in  obtaining  such  results  can  reasonably  be 
expected  by  rephrasing  the  problem  in  the  generalized  density  operator 
formalism  and  making  better  use  of  the  available  tools  of  modern 
functional  analysis.  Preliminary  studies  along  this  line  are  presently 
being  conducted  here  at  Stanford  Research  Institute  (ISR- 18453  1  -  1 10). 
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I 


III.  THE  T-F  AND  THE  T-F-D  STATISTICAL  MODELS  FOR  MONO- 
ATOMIC  CRYSTALS3 

A .  Introduction 

The  semiclassical  statistical  approach,  introduced  by 
Thomas  and  Fermi  in  1927  and  developed  since  by  many  others,  is 
very  simple  conceptually  but  its  range  of  validity  and  proper  estimates 
of  the  errors  involved  (as  compared  to  the  Hartree,  Hartree -Fock, 
and/or  the  exact  solution  for  a  given  Hamiltonian)  have  not  yet  been 
established,  although  significant  progresses  have  recently  been  made 
in  answering  such  questions,  in  particular  by  March,  Kirzhnits  and 
others.  A  discussion  of  the  foundations  of  the  statistical  method  is  best 
conducted  within  the  frame  of  the  generalized  density  operator  theory 
and  will  be  given  elsewhere.  Here  we  shall  give  only  an  heuristic  and 
a  variational  derivation  of  the  basic  equation  of  the  T-F  and  T-F-D 
models  and  shall  compare  the  resulting  equations  of  state  for  solids 
with  those  already  obtained. 


B .  Description  of  the  Models  and  Heuristic  Derivation  of  the 

Basic  Relation 

In  the  simplest  statistical  models  a  solid,  represented  by 
a  system  of  electrons  in  a  positive  lattice  background,  is  treated  in  the 
following  manner.  The  kinetic  energy  of  the  nuclei  is  neglected,  and 
the  validity  of  the  static  approximation  is  assumed;  furthermore,  one 
argues  heuristically  that: 

1)  The  electron  system  can  be  treated  locally  as  a  free 
electron  gas--with  exchange  in  the  T-F-D  model-- 
raised  at  some  position  dependent  electrostatic 
potential  V(r). 
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2)  The  electron  system  is  in  a  bound  state  and  in  statis¬ 
tical  equilibrium,  i.  e.  ,  no  energy  can  be  gained  by 
transferring  an  electron  from  one  point  to  another  in 
the  solid. 

3)  The  electrostatic  potential  of  the  electrons  can  be 

taken  as  that  of  an  equivalent  (charge  wise)  continuous 

charge  distribution  -ep  (*)  satisfying  the  classical 

Poisson  equation  with  the  usual  continuity  conditions 
+  dV 

for  V(r)  and  €  - - - -  across  any  boundary 

n  dn 

surface . 

AVe(r)  =  4ltepe(r) 

(1) 

+  b  .  c  . 

4)  The  total  energy  of  the  electron-nuclei  system  can 

be  obtained,  in  the  T-F  model,  as  the  sum  of  the  free 
electron  kinetic  energy  E^  and  the  "classical"  elec¬ 
trostatic  potential  energy  of  the  system  Ec|  .  In  the 
T-F-D  model  the  Pauli  correlation  effects  are 
approximately  taken  into  account  by  adding  the  free 
electron  gas  exchange  energy  E^ 

From  these  assumptions,  one  derives  the  basic  relation  of  the 

e  ^ 

model,  namely,  the  relation  between  p  (r)  and  the  total  electrostatic 
potential  V(r) 

V(?)  =  V  (?)  +  V  (?)  +  V  J?) 

n  e  ext 

where  V  is  the  potential  due  to  the  nuclei,  V  that  due  to  the  eier- 
n  e 

trons  and  V  '  any  external  one. 
ext 


❖  V  is  assumed  null  in  the  sequel, 
ext 


157 


From  1)  and  2)  it  follows  that  the  energy  of  an  electron  at  a  point  r 
is  given  by 


2  -» 
P  (r) 

2m 


-  eV(r)  -  —  p(r) 


Since  one  is  considering  bound  states,  one  must  also  have  the  relation 


-£r  -  eV'r)  •  w  p<r> s  -  eV0 


(2) 


where  -eVQ,  the  maximum  potential  energy  of  an  electron  at  r,  is 
position  independent  (from  condition  c)  ). 


The  maximum  momentum  p  (r)  of  a  free  electron  gas  is  re¬ 
max 


lated  to  the  electron  density  by 


-»  2  1/3  e 

Pmax'r>  *  <3"  >  ^ 


1/3 


Combining  these  two  relations,  one  obtains  the  basic  relation  of  the 
T-F  and  T-F-D  models: 

1  2  -> 

-j-  P  (r)  -  P  (r)  =  e(V -  V  ) 

2m  max  fffi  max  o 


or 


2  2/3 


2  2  /3  ft  e 

=(V-V  1  <3ff  >  -ZST  p 


_3_ 

V 


1/3  3  e'/J 

e  p 


(3) 


5  e2/3  4  e 

T  V  "  T  V 


1/3 


or 


with 


_  (2em) 

A  =  - - 


pC(r)  -  A  ( 
3/2  I 


2  - 
a  +  [  a  +  ( V  - V  )J  z 


2  3 
3  ff  ft 


Je 


3/2 


(4) 


5a. 


a  =  (2em)2  ~  - 


2a 


(15eak) 


(4') 


(For  the  T-F  model,  one  sets  a  -  0.) 
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This  relation,  of  course,  holds  only  in  regions  where 


2 

a  +  (V-V  )  *  0 
o 

and  in  regions  of  zero  electron  density,  the  Poisson  equation  is  replaced 

e  -+ 

by  Laplace's  equation.  Knowing  p  (r)  and  V(r),  the  total  energy  of 
the  system  is  obtained  from: 


E 


I  P  (r)  f  k(r)dr  + 


/ 

J 


e 

P  (r) 


C  (?)dr 


+ 


+  V 

e 


e  -> 
p  (r) dr  +  E 


n-n 


(5) 


C.  Variational  Derivation  of  the  Basic  Relation 


It  is  of  interest  to  note  that  the  basic  relations  (3)  or  (4)  can 
also  be  obtained  from  a  variational  principle,  namely  by  minimizing  the 

total  energy  E--as  given  by  (5)--under  variation  of  the  electron  density 

e  ^ 

p  (r)  for  fixed  boundaries  {and  so  fixed  nuclei)  and  subject  to  the  sub¬ 
sidiary  condition 

J  pG(r)dr^  =  /T  =  const  (6) 

<>r .  fixed  total  number  of  electrons.) 


is 


Introducing  a  Lagrange  multiplier  V  ,  the  extremum  condition  is 


6  [  E  +  e  V  ]  -  0  . 


(7) 


Using  the  expressions  given  previously  for  €  and  €  and  the  relation 

K  X. 


V  (r) 
e 


--  [Jm** 

J 


pe(h^ 

(i.e.  ,  Poisson's  equation)  one  can  write  the  energy  in  the  form: 

.5/3 


(8) 


=  /[<*, 


e  ~ e,->.4  /31  <♦  e  ff 

p  (r)  -  O^p  (r)  ]  dr  +  y  J  J 


lpe(r)pe(r') 


drdit 


I 


r’  -  r 


I 


-  e 


f  e  “► 

/  V  (r)p  ( r)dr  +  E 
J  n  n-n 


(9) 
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and  Eq  (7)  gives 


1 1“ 


5  e 

k  T  P 


2/3 


(r) 


4  e  ->  4  f  4, 

-  —  axP  (r)  ]  dr  -  e  J  (vn“v0)dr 


or 


L  U/H l  d?d?.  +  4-  \\£&-  d?d? '=  o 
:  j  J  if .  ?i  2  JJ  i?. .  ti 


i 


5  e2/3  4  e1/3 

T  V  -  T  V 


-  e( V - V  )>  dr  =  0 

o  ‘ 


In  view  of  the  bound  state  condition  (2),  it  follows  that  the  Lagrange 
multiplier,  V  ,  is  the  maximum  electron  potential  and  that 


e(V " Vq)  =  T  V 


.2/3  4 


1/3 


,  “  P 
3  x 


Q.E.D. 


D.  Some  General  Remarks  on  the  T-F  and  T-F-D  Models 

The  basic  quantitites  in  the  statistical  models  are  the  elec- 
tron  density  p  (r)  and  the  electrostatic  potential  V(r)  =  V^(r)  +  VMr) 
+  eventually  an  external  potential.  An  heuristic  argument  or  an  equiva¬ 
lent  variational  procedure  leads  to  the  basic  nonlinear  relations  Eqs  (3) 
or  (4)  of  the  T-F  and  T-F-D  theories;  the  assumption  of  the  validity 
of  Poisson's  equation  and  the  usual  electrostatic  boundary  conditions 
complete  the  description  of  the  models. 

In  discussing  this  approach,  one  can  distinguish  two  groups  of 
specific  assumptions,  besides  those  leading  to  the  corresponding 
Hamiltonian  of  the  proper  quantum  mechanical  approach: 

1)  The  exact  energy  expression  for  the  given  Hamiltonian 
is  replaced  by  a  simpler  one  where  not  only  the  cor¬ 
relation  terms  are  omitted,  but  also  the  kinetic  energy 
is  simplified  by  neglecting  the  effect  of  density  gradi¬ 
ents  (off  diagonal  part  of  Tfx'j.Xj),  the  1st  order 
generalized  density  matrix). 
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2)  The  description  of  the  system  in  terms  of  the  wave 
function  *(xr...)  or  the  reduced  density  matrices 

“(xl  >  xj)>  ^(xix2^  *s  reP^aced  hy  a  description  in 
terms  of  p(T)  alone,  i.e.,  essentially  in  terms  of  the 
diagonal  part  only  of  the  first  order  density  matrix, 

r  (x^ ,  x  l). 

C 

3)  The  electron  density  p  (r)  is  calculated  on  the  basis 

of  heuristically  derived,  essentially  classical  equations, 
rather  than  by  solving  the  proper  Schrodinger  equation. 

The  difficulties  in  assessing  the  exact  significance  of  this  approach 
lie  essentially  in  the  last  two  assumptions.  A  proper  foundation  of  the 
statistical  method  would  require  not  only  an  estimate  of  the  neglected 
term's  in  the  energy  expression,  but  also  a  derivation  of  Eqs  (1),  (3),  or 
(4)  as  a  meaningful  approximation  to  the  quantum  mechanical  ones. 
Although  several  authors  have  recently  made  significant  progress  in 
this  direction,  no  complete  solution  of  the  problem  is  yet  available.  It 
should  also  be  noted  that  the  statistical  approach  (T-F,  T-F-D  suitably 
modified)  is  not  necessarily  limited  to  the  accuracy  of  the  Hartree-Fock 
approximation  but  could  very  well  offer  in  specific  problems,  a  better 
approximation  to  the  exact  solution,  especially  in  view  of  the  difficulties 
in  obtaining  any  degree  of  self-consistency  in  H-F  calculations  for 
solids . 

E .  The  T-F  and  T-F-D  Models  for  Simple  Monoatomic 

Crystals 

1.  General  Equations 

For  a  crystal  of  lattice  given  by 

R,  .  =  Ln.a,  r..  =±,0  integers 

(nl  ill  i 
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and  such  that  all  lattice  sites  are  occupied  by  identical  atoms  of  atomic 
number  Z,  the  energy  per  atom  in  the  T-F-D  method  can  be  written 
as: 


E 


n 


lim 
N— > 00 


e  ->  •+  ■> 

P  (rKk(r)dr  + 


/ 

N  0 


e  -►  -> 

P  (r)Cx(r)dr 


r  e  2  2 

y  I  Zep  (r)  i  y  y1  Z  e _ 

(m)  " ?l  ^ (m)(n>  !?(m)  -  *(n) ' 

+  t  ff  Pe(?)Pe(r)  d*d?l 

NO  I -  r  |  J 

_  ^  ^ 

where  O  is  the  primitive  cell  of  volume  =  (a  ,  a  ,  a  ).  (For  the  T-F 

A  Lt  J 

expression  one  should  take  -  0.  ) 

Using  the  translational  invariance  of  the  crystal 

p  (r  +  R<m)>  =  P  (r)  V(r+  R(n)r  V(r) 


to  simplify  this  energy  expression  one  obtains: 


and  the  boundary  value  problem  giving  the  electron  density  reduces  to  a 
boundary  value  problem  in  a  single  cell  Cl  . 

e 

For  the  simple  monoatomic  solid  considered  here  p  (r)  is  to  be 
obtained  as  a  solution  of  the  system  of  equations: 
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(11) 


r  5  eL '  ^  4  e-*l/3 

e[V(r)  -  Vq]  =  —  «kP  (r)  -  —  O^p  (r)  7 


A[  V(r)  -  Vq]  =  41T  ep  (r)  , 


with  the  boundary  conditions 


lim  rV(r)  =  lim  r[V-V  1  =  eZ 

x*  p  q  r  — » o  L  oJ 


(12) 


=  0 


Sq  =  surface  of  the  cell 


The  first  boundary  condition  states  that  the  total  potential  converges  to 
that  of  the  nuclei  +  Ze  as  r — ^o  ;  the  second  follows  from  the  lattice 
translational  invariance,  the  continuity  requirement  on  and  the 
inversion  symmetry  of  0  and  V.  One  should  note  that  the  translational 
invariance  automatically  insures  the  continuity  of  V(r)  across  S  (if  V 
is  bounded  on  S). 

2  .  The  Sphere  Approximation 

The  "sphere  approximation"  consists: 

a)  in  replacing  the  W-S  polyhedron  cell  0  by  a  sphere 

of  equal  volume,  the  so-called  Wigner-Seitz  sphere  of  radius  r  given 

s 


=  f^i-  til1 

[  4ff  o 


3vaU/3 


(13) 


where  v  =  atomic  volume  ; 
a 


b)  in  assuming  that  the  electron  density  and  elec¬ 
trostatic  potential  are  spherically  symmetric 


Pe  =  pG(r)  and  V  =  V(r) 


r  =  I  r  I  ;  and 
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c)  in  neglecting  the  overlap  of  the  "equivalent 
spheres"  in  calculating  the  energy  of  the  system. 

G 

With  these  supplementary  assumptions  the  equations  for  p  take 
the  form 

e[V(r)  -  VJ  =  ~  <*kpe(r)2/3  -  ~  cy>e(r) 1/3 


1  d  e 

- - y  (  rV(r)  )  =  4lTp  (r) 

r  dr 


with  the  boundary  conditions 


(14) 


lim  ,  ,  .  , 

r— *o  r[  V(r)  -  V  ]  =  eZ 


dV 

dr 


=  0 


(and  V(r  )  =  0  ) 
s 


r=r 


(15) 


and  the  energy  expression  reduces  to 

r 

,5/3 


En  *  4n 


J 


-axpe 


s  e, 


4/3 


2  . 

r  dr 


.  4  ne°f  5.-2  V  [ 

Jr  J  |r'  -  r| 


(16) 


e  2 
p  (r)r  dr, 


i.  e.  ,  to  that  of  an  "equivalent  spherical  atom.  "  The  only  difference  be¬ 
tween  a  free  atom  problem  and  that  of  a  solid  in  the  sphere  approximation 

g 

is  that  in  the  later  case  p  (r^)  is  not  equal  to  zero  in  general,  so  the 
pressure  at  the  atom  boundary  due  to  the  electrons  is  also  different 
from  zero  (compressed  atoms). 


Note  tha  t  the  overlap  correction  terms  neglected  are  of  the  form 

r 


2  2 
e  Z 


<m)  2  I  R(m)| 


-  477  eZ 


/ 


S  pe(r)r2dr 

I  a  -t 

I  (m) 


+ 
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(16) 


+  (4ff ) 


II 


S  ^  I  \  ^  /  I  \ 

r.2r2  P  (r)P  (r  ) 


r1  -  r  +  R 


dr 'dr 


(m) 


and  were  included  in  the  calculation  of  the  energy  of  the  uniform  elec¬ 
tron  gas  in  a  positive  lattice  background. 

F.  Some  Remarks  on  the  Sphere  Approximation  and  on  the 
Inclusion  of  Correlation  Effects  in  the  T-F-D  Model 

Clearly  the  sphere  approximation  is  in  no  way  justified  for 

solids  and  is  made  solely  for  convenience  to  reduce  the  3 -dimensional 

c 

boundary  value  problem  giving  p  to  a  1 -dimensional  one.  Some  of 
the  physical  consequences  of  this  approximation  are  that: 

1)  The  crystallographic  structure  of  the  solid  is  almost  com¬ 
pletely  washed  away,  the  volume  of  the  cell  being  the  only  remaining 
geometrical  parameter. 

2)  Since  the  solid  is  replaced  by  a  collection  of  neutral  spheres 
with  spherical  symmetric  charge  distributions,  no  energy  is  involved 

in  forming  the  solid- -only  in  compressing  the  individual  atoms,  and 
this  model  cannot  offer  any  resistance  to  shear  stresses. 

3)  The  extension  of  such  a  model  to  nonsimple  crystals,  in 
particular  to  polyatomic  ones,  can  be  done  only  in  a  very  artificial  way 
(such  as  the  smoothing  technique)  or  by  introducing  "effective  mono- 
atomic  models,  "  i.e.  ,  through  largely  artificial  and  arbitrary  averag¬ 
ing  procedures . 

4)  Calculation  of  the  electronic  structure  of  the  solid  in  this 
approximation  is  almost  meaningless. 

If,  at  the  time  of  the  development  of  the  T-F,  T-F-D  models, 
the  available  mathematical  tools  seemed  to  leave  little  hope  of  dealing 
successfully  with  the  3 -dimensional  problem,  the  recent  developments 
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in  mathematical  techniques  for  handling  such  nonlinear  problems- -in 
particular  variational  ones --have  completely  altered  the  picture,  and 
it  appears  that  approximate  solutions  and  proper  error  estimates  could 
be  obtained  for  the  actual  problem. 

This  is  one  direction  along  which  further  research  at  this  time 
is  possible.  Within  the  limits  of  the  statistical  approach  it  would  allow 
one  to  obtain  a  more  complete  and  realistic  model  for  monoatomic  and 
polyatomic  solids  at  high  pressures  including  the  pressure  dependence 
of  the  elastic  coefficients,  dielectric  tensor  and  degree  of  anisotropy; 
possible  high  pressure  phase  transformations  could  be  better  handled 
and  eventually  the  electronic  structure  of  different  types  of  materials 
could  also  be  studied  over  an  extended  density  range  by  solving  the 
Bloch  problem  with  the  3 -dimensional  statistical  crystal  potential. 

Correlation  effects  can  be  introduced  in  the  T-F-D  model  in  the 
same  approximate  way  as  the  exchange  effect  in  the  T-F  model,  namely, 
by  adding  to  the  energy  expression  a  term  arising  from  the  correlation 
energy  of  the  locally  free  electron  gas: 

e  [pe*)7  [Pe]dt. 

corr  J  corr  1 

This  has  been  done  by  Gombas,  Lewis,  and  Erma.  The  resulting  basic 

0 

relation  between  p  and  V(r)  becomes  quite  involved,  and  a  straight¬ 
forward  numerical  treatment  in  the  line  of  those  applied  to  the  T-F-D 
model  is  almost  prohibitive.  Also  the  value  of  this  extension  of  the 
statistical  method  is  very  doubtful  for  it  still  leaves  out  completely 
the  inhomogeneity  corrections*  which  can  be  at  least  of  the  same 
order  as  the  Coulomb  correlation  one. 

A  second  direction  along  which  profitable  work  on  the  statistical 
approach  could  be  done  presently  is  in  obtaining  meaningful  extensions 

*  These  are  essentially  kinetic  energy  correction  terms  that  take 

partly  into  account  the  effect  of  potential  (or  electron  density)  gradients . 
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of  the  T-F-D  models;  these  should  incorporate  both  inhomogeneity 
corrections  and  correlation  effects  in  a  consistent  fashion--at  least  to 

the  same  order  in  a  perturbation  expansion- -and  should  include  a 

10 

practical  computational  scheme.  Various  formal  extensions  have  been 
proposed  in  the  last  few  years  by  several  authors,  but  little  numerical 
work  is  available  outside  free  atom  calculations. 

These  different  schemes  are  being  presently  evaluated,  partly 
in  conjunction  with  basic  work  on  the  statistical  and  generalized  density 
operator  theory,  and  we  shall  eventually  report  later  on  their  possible 
use  in  deriving  improved  equations  of  state  for  solids. 
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